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Abstract 



We derive the imaginary part of the potential NRQCD Hamiltonian up to 
order 1/m 4 , when the typical momentum transfer between the heavy quarks is 
of the order of Aqcd or greater, and the binding energy E much smaller than 
Aqcd- We use this result to calculate the inclusive decay widths into light 
hadrons, photons and lepton pairs, up to 0(mv 3 x (AQ CD /m 2 , E/m)) and 
0(mv 5 ) times a short-distance coefficient, for S- and P-wave heavy quarko- 
nium states, respectively. We achieve a large reduction in the number of 
unknown non-perturbative parameters and, therefore, we obtain new model- 
independent QCD predictions. All the NRQCD matrix elements relevant to 
that order are expressed in terms of the wave functions at the origin and six 
universal non-perturbative parameters. The wave-function dependence factor- 
izes and drops out in the ratio of hadronic and electromagnetic decay widths. 
The universal non-perturbative parameters are expressed in terms of gluonic 
field-strength correlators, which may be fixed by experimental data or, alter- 
natively, by lattice simulations. Our expressions are expected to hold for most 
of the charmonium and bottomonium states below threshold. The calculations 
and methodology are explained in detail so that the evaluation of higher or- 
der NRQCD matrix elements in this framework should be straightforward. 
An example is provided. 
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I. INTRODUCTION 



Heavy Quarkonium is characterized by the small relative velocity v of the heavy quarks in 
their centre-of-mass frame. This small parameter produces a hierarchy of widely separated 
scales once multiplied by the mass m of the heavy particle: m (hard), mv (soft), mv 2 
(ultrasoft), .... In general, we have E ~ mv 2 <C p ~ mv <C m, where E is the binding 
energy and p the relative three-momentum. 

The use of NRQCD [1] allowed a factorization of the physics due to the scale m from the 
one due to smaller scales. Moreover, it allowed the description of heavy quarkonium inclusive 
decays into light fermions, photons and leptons, in terms of matrix elements of local 4-quark 
operators, in a systematic way. These 4-quark operators are of two types: colour-singlet and 
colour-octet operators. The matrix elements of the colour-singlet operators can be related in 
a rigorous way with quantum- field theory defined quarkonium wave functions [1]. Intuitively, 
these wave functions should be related with the wave functions that appear in a Schrodinger- 
like formulation of the bound-state system, namely two heavy quarks interacting through 
a potential. On the other hand, the colour-octet ones were thought to have no parallel 
in such formulation. In either case, even though there had been a lot of relevant work in 
obtaining the QCD potential in terms of Wilson loops [2], it was not known how to obtain 
the systematic connection between NRQCD and a Schrodinger-like formulation in the non- 
perturbative case, or even whether it existed and, if so, under which circumstances. Even 
in the perturbative case, for which expressions for the potential existed at lower orders in 
the past [3] , a clean and simple derivation of such Schrodinger-like formulation incorporating 
perturbative ultrasoft gluons was not clear once higher order calculations in a s were required. 

The observation that NRQCD still contains dynamical scales that are not relevant to the 
kinematical situation of the lower-lying states in heavy quarkonium (those with energy scales 
larger than the ultrasoft scale) [4] (see also [5]), paved the way towards the resolution of the 
questions above. Indeed, it was realized that further simplifications occur if we integrate 
them out, and the resulting effective field theory was called pNRQCD [4]. The degrees of 
freedom of pNRQCD depend on the interplay between the characteristic scales of the given 
non-relativistic system, namely E, p and the momentum transfer k, and the characteristic 
scale of non-perturbative physics in QCD, which will be denoted by Aq CD . Therefore, how a 
Schrodinger-like formulation develops, and thus how the NRQCD 4-fermion matrix elements 
will show up within this framework, depends on the specific kinematic situation considered. 

When the typical momentum transfer k is much larger than Aqcd, k ~ p ^> E > Aqcd, 
the pNRQCD Lagrangian [4,6] contains not only the singlet field, which is also present in 
the Schrodinger-like formulation, but also the octet field, ultrasoft gluons and light quarks. 
The matching from NRQCD to pNRQCD (integration of the soft scale) can be done in 
perturbation theory. In nature, this situation is relevant to the Y(IS') and the t-i production 
near threshold. If in addition E ^> Aqcd, we are entirely in the weak-coupling regime 
(E ~ ma 2 , p ~ k ~ ma s ) where non-perturbative effects can be parameterized by local 
condensates [7]. In this regime pNRQCD has been used to obtain the complete set of 
logarithmic corrections to the QCD static potential at three loops [8], the complete set of 
logarithmic corrections to the very heavy quarkonium spectrum at 0(m« s 5 ) [9] (see also 
[10]), the resummation of logs at the same order [11,12] and, very recently, the (almost) 
complete spectrum of very heavy quarkonium at 0(ma s 5 ) [13]. We can still use the same 
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pNRQCD Lagrangian for systems with E > Aqcd- Then, however, some of the calculations 
in pNRQCD cannot be carried out perturbatively and the non-perturbative effects can no 
longer be parameterized by local condensates (see [7,6]). 

When the typical momentum transfer k > Aqcd and the binding energy is small, namely 
E <C Aqcd, the degrees of freedom of pNRQCD are the singlet field and pseudo-Goldstone 
bosons (pions), if hybrids and other degrees of freedom associated with heavy-light meson 
pair threshold production develop a mass gap of C(Aq CD ), as it is assumed in Refs. [14,15,6] 
and in what follows. If we ignore Goldstone bosons, which play a negligible role in the 
present analysis, we recover the celebrated Schrodinger-like picture of quark and antiquark 
interacting through a potential. Therefore, the pNRQCD Lagrangian reads [14,6]: 



where h is the pNRQCD Hamiltonian, to be determined by matching to NRQCD. In general, 
one should be able to obtain the binding energies and the total decay widths from the real 
and imaginary parts of the complex poles of the propagator. At the accuracy we are aiming 
at in this paper the total decay width of the singlet heavy-quarkonium state may be defined 
as: 



where \n, I, s,j) are the eigenstates of the real part of the Hamiltonian h. 

In this paper we will be concerned with this situation and will consider in full detail not 
only the calculation in the general case (A) Aqcd < k (Section III), but also the particular 
situation (B) Aq C d "C k (Section V): 

A) Aqcd is smaller than or of the order of k. In this case, the (non-perturbative) matching 
to pNRQCD has to be done in a single step. This case has been developed in a systematic 
way in Refs. [14,15]. As a consequence the complete set of potentials up to order 1/m 2 could 
be finally calculated [14,15], including a 1/m potential, which had been missed so far and 
completing (and in some cases correcting) the previous expressions obtained in the literature 
[2] for the 1/m 2 potential. Most of the charmonium and bottomonium states below threshold 
are expected to be in this situation. 

B) Aqcd is much smaller than the typical momentum transfer k. In this case, the degrees 
of freedom with energy larger than or similar to k can still be integrated out perturbatively. 
This leads to an intermediate EFT that contains, besides the singlet, also octet fields and 
"ultrasoft" gluons (meaning gluons with energies < Aqcd here) as dynamical degrees of 
freedom [4,6]; it has the same Lagrangian as pNRQCD in the weak coupling regime. We will 
call this EFT pNRQCD'. 1 The octet and "ultrasoft" gluon fields are eventually integrated 
out by the (non-perturbative) matching to pNRQCD [6]. 

In either case, it remained to be seen how the matrix elements of the 4-fermion operators 
are encoded in this formulation. This was especially needed for the octet ones since, as 
mentioned before, it was thought that they could not be accommodated in a Schrodinger-like 



Note the change of name with respect to section 5 of [6]. 
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formulation. However, in [16], we have shown that, by using pNRQCD, it is indeed possible 
to relate the matrix elements of the colour-octet operator with the wave function at the 
origin and additional bound-state independent non-perturbative parameters. This was done 
for the specific case of P-wave quarkonium decays. Here, we will apply the same method to 
express all the NRQCD matrix elements relevant to inclusive S'-wave quarkonium decays into 
light hadrons, photons and lepton pairs at 0(c(a s (m))mv 3 x (AQ CD /m 2 , E/m)) (c(a s {m)) 
being a function of a s (m) computable within perturbation theory). This reduces the number 
of unknown parameters for the total decay widths of charmonium and bottomonium states 
below threshold by roughly a factor of 2, which allows us, in turn, to formulate several new 
model-independent predictions. Particularly important is the fact that our formalism allows 
the physics due to the solution of the Schrodinger equation, which appears entirely in the 
wave-function, to be disentangled not only from the short-distance physics at scales of 0(m), 
but also from the gluonic excitations with an energy of 0( Aqcd) • As a consequence, the wave- 
function dependence drops out in the ratio of hadronic and electromagnetic decay widths. 
For this class of observables the reduction in the number of non-perturbative parameters in 
going from NRQCD to pNRQCD is even more dramatic, since only the (six) non-perturbative 
universal parameters appearing at this order in pNRQCD are needed. 

Finally, we would like to mention the dynamical situation when the binding energy is 
positive and of the same order of magnitude as the momentum transfer k, namely when 
E > Aqcd ~ k. In this case degrees of freedom with energy ~ Aqcd cannot be integrated 
out. States close to and beyond heavy-light meson pair threshold are expected to be in this 
situation. The results of this paper do not apply, in principle, to this case. 

The paper is arranged as follows. Section II reviews some aspects of NRQCD that 
are relevant to the rest of the paper. Section III provides a detailed description of the 
computation of the "spectrum" of NRQCD, in particular the ground state, in the 1/m 
expansion in the general case. It is meant for the reader interested in learning the techniques 
involved in this type of computations. The description of pNRQCD, its power counting 
and the relation between the computation of Section III and the Hamiltonian in pNRQCD 
are given in Section IV. Section V provides a detailed description of the matching between 
pNRQCD and NRQCD in the particular case E <C Aq C d <^ k. This section may help 
the reader who is not willing to go through the general case in Section III, but still wants 
to get a flavour of the kind of calculations we are performing. Section VI summarizes 
our results. The reader who is only interested in our final results and wants to skip any 
computational detail may jump directly to this section. Section VII displays some model- 
independent predictions that follow from our results. We finally draw our conclusions in 
Section VIII. A number of appendices complement the main body of the paper. Appendix 
A recalls the 4-fermion NRQCD operators at (9 (1/m 4 ). Appendix B gives the general formula 
relating an arbitrary NRQCD matrix element with the computation in pNRQCD. Appendix 
C gives the leading log renormalization group running of the imaginary parts of the 4-fermion 
NRQCD operators matching coefficients. Appendix D shows how to deal with ill-defined 
products of distributions within dimensional regularization. Appendix E shows how unitary 
transformations can relate different forms of the pNRQCD Hamiltonian. 
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II. NRQCD 

NRQCD is obtained from QCD by integrating out the heavy-quark mass scale m [1]. The 
NRQCD Lagrangian can be written as follows: 

£nRQCD = £ g + Aight + A-f + (3) 

where £ g involves only gluon fields, £n g ht involves light-quark and gluon fields, and An-f 
are the terms in the Lagrangian with 2n heavy-quark fields. 

The NRQCD Lagrangian can be organized as a series expansion in a s (m) and in the in- 
verse of the heavy-quark mass 1/m. Powers of a s (m) are encoded into the Wilson coefficients 
of NRQCD. 

In this paper, we aim at a description of heavy quarkonium inclusive decays into light 
hadrons and electromagnetic decays, whose appearance is due to the imaginary terms of 
the NRQCD Lagrangian. It is convenient, then, to split the Lagrangian into the Hermitian 
(real) and the anti-Hermitian (imaginary) part: 

£nrqcd = R- e ^nrqcd + i Im £nrqcd, (4) 

where 

Re £nrqcd = £g + flight + A-f + Re A-f , (5) 

and 

Im£ NRQCD = Im£ 4 _ f . (6) 

The operators responsible for heavy-quarkonium decays are the NRQCD 4-fermion operators 
whose matching coefficients carry an imaginary part. For our purposes, it is sufficient to 
consider either dimension 6 or dimension 8 4-fermion operators: 

Im £ NRQCD = Im £ 4 _ f = Im C£$ + Im £™ d = 6 + Im Cfzl + Im £™ d = 8 . (7) 

With the superscript EM, we indicate operators responsible for the electromagnetic decays. 
More explicitly, we have 

Im £ « = MM 0l( . Sl) + 0l( 3 Sl) + Os cs o) 

m 2 m 2 m 2 



Im/ 8 ( 3 ^) 



m 



2 



Os( i S 1 ), (8) 



Im£ EM^ 6 = Wem(%) 0em(1 ^ o) + Wem^ (35 } (9) 
m 2 m 2 



"I" | - ' 1 V ' - ! I ' H "U.l j 

m 4 m 4 m 4 



+ Img 1 (ZS 1 *D 1 ) ^ + ^ ^ ^ ^ ^ 



Ttyi r EMd=s Im /em! 1 -Pi) n ,i pU ImW^o) n , 3 p \ , Im /EM( 3 -Pi) ^ / 3p n 

im£. 4 _ f — CemI, -mJ H 2 ^EMV -H)J H 3 V-emK r\) 

m 4 m 4 m 4 

+ Im/ E M( 3 P 2 ) 0EM 3 P2) + ImW^o) + Im,EM( 3 ^) pEM , 

m 4 m 4 m 4 

+ lm9EM{3S y Dl) V EM ( 3 Su 3 D 1 ). (11) 
m 4 

The definitions of the hadronic operators can be found in [1]. For ease of reference, we recall 
them in Appendix A, where we also give the definitions of the electromagnetic operators. 

The distinction between hadronic and electromagnetic operators is somewhat artificial. 
In general the 4-fermion operators listed in Eqs. (Al)-(A18) are all the dimension 6 and 8 
operators needed to describe decays into light hadrons and/or hard electromagnetic parti- 
cles. The information needed in order to describe decays into hard electromagnetic particles 
is encoded into the electromagnetic contributions to the matching coefficients. The electro- 
magnetic operators defined in [1] arise from singling out the operators accompanying the 
matching coefficients whose imaginary parts correspond to pure electromagnetic decays and 
inserting into them the QCD vacuum (|vac)(vac|). This insertion guarantees that when cal- 
culating with these operators in NRQCD, no contamination from soft strong interactions will 
occur. Hence, the electromagnetic operators encode all the relevant information needed in 
order to calculate the quarkonium total decay width to electromagnetic particles only. How- 
ever, one might also be interested in the decays to hard electromagnetic particles and soft 
light hadrons. In this case, the complementary to the above projector, namely 1 — |vac)(vac| 
should be considered. In this paper, however, we will restrict to the processes, and therefore 
to the operators, originally considered in [1]. 

The Hermitian piece of the NRQCD Lagrangian can also be written in a 1/m expansion: 
Re£ = + + ^Re£( 2 ) + • • • . (12) 

At order 1/m the different pieces of Eq. (5) read 

C,_, = ,^D + 5! + c F9 ^y + X *\iD - 5! - c F9 ^} X , (13) 



flight = ^qjipqj, 
j'=i 

Re £ 4 _f = 0, 

where ip is the Pauli spinor field that annihilates the fermion and \ is the Pauli spinor field 
that creates the antifermion, iD° = id - gA°, iD = iV + gA, B* = -e ijk G jk /2; for later 
use, we also define E* = G j and [D-, E] = D ■ E — E • D. The chromomagnetic matching 
coefficient cp is known at next-to-leading order and its value can be found in [17]. Concerning 
the explicit expression of the 0(l/m 2 ) Lagrangian, see Ref. [15] for the operators without 
light quarks and Ref. [18] for the operators including light fermions. 
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III. THE NRQCD "SPECTRUM" IN THE 1/M EXPANSION 



We assume to be in the situation Aqcd < tuv in which the matching to pNRQCD cannot 
be performed within a perturbative expansion in a s . Nevertheless, it can be done by assuming 
an expansion in 1/m, within the Hamiltonian formalism of [14,15], to which we refer for 
further details. We may divide the procedure into three steps: 

1) The spectrum of the NRQCD Hamiltonian, made of quarkonium and gluonic excita- 
tions between heavy quarks, is evaluated order by order in 1/m starting from the static 
configuration. This will be done in Sees. Ill A-III E. 

2) The quantum-mechanical matrix elements are expressed in terms of gluonic field cor- 
relators. This will be done in Sec. Ill F. 

3) The excitations of order mv 2 are identified as the degrees of freedom of pNRQCD. The 
matching to pNRQCD is performed by integrating out the excitations of order Aqcd 
and mv. This will be done and discussed in Sec. IV. 



A. The NRQCD Hamiltonian 

The NRQCD Hamiltonian without light fermions has been worked out up to 0(l/m) in 
Ref. [14] and up to 0(l/m 2 ) in Ref. [15], to which we refer for the explicit expressions. In 
the following we will consider the inclusion of light fermions. 

The inclusion of light fermions produces new terms in the Hamiltonian of pure gluody- 
namics. In the static limit, we have: 



n/ r 

H (o) = #(o)( n/ = 0)-J2 iD -yqj. 

.7=1 J 



(14) 



The next corrections in the Hamiltonian, due to light fermions, appear at (9 (1/m 2 ) and have 
been considered in Ref. [18]. We will not need their explicit expressions in this paper. We 
will only need the expressions of the Hermitian part of the NRQCD Hamiltonian up to order 
1/m: 

ReH= — = -— J rfW (B 2 + gc F CT-B)^ + — J d 3 x X f (D 2 + g c F <r • B) X , (15) 
and of the imaginary part of the NRQCD Hamiltonian up to order 1/m 4 : 

lmH = — + — , (16) 

m A m 4 

where Im#( 2 ) = Imtf^, ImtfW = lmH^} f and 

-^ = -jd^(£f_} + £™n, (17) 

n 4-f _ 



m 4 



J rf 3 x (jCfZ 8 f + Cf M /= 8 ) . (18) 
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The Gauss law, constraining the physical states |phys), reads: 

D • IT |phys> = g{^T a ^ + x ] T a x + g7°T a g)| P hys>, (19) 

where II a is the canonical momentum conjugated to A a . In Ref. [19], general details about 
Hamiltonian quantization can be found and in Refs. [14,15] details specific to our case. 

B. The NRQCD spectrum at C(l/m 3 ) 
Let us call H = + Hj the NRQCD Hamiltonian, if (°) being its static part and 

Hj = ++.... 20 

m m z 

We call |n;xi,x 2 )(°) the eigenstates of H , the corresponding eigenvalues, |n;xi,x 2 ) the 
eigenstates of H, and E n the corresponding eigenvalues within a strict expansion in 1/m. 
This means that they satisfy the analogue of the Schrodinger equation: 

H\ m xi, x 2 ) = J dV^V^n; x' l5 x'^K, x' 2 ; p^, p' 2 )5^ (xi - Xl )^ 3 )(x 2 - x 2 ). (21) 

With n we indicate a generic set of conserved quantum numbers. Note that the heavy- 
quark positions x x and x 2 are conserved quantities only with respect to the zeroth order 
Hamiltonian H . The states |n;xi,x 2 ) are normalized according to 

(m; xi, x 2 |n; yi , y 2 ) = 5 nm 5 (3) ( Xl - yi)5 (3) (x 2 - y 2 ), (22) 

and we define 

^ /2 (Yi, y 2 ; Pi, p 2 )5 (3) (yi - Xi)6®(y 2 - x 2 ) = <°> (n; yi , y 2 |n; x 1; x 2 ). (23) 

The above three equations (21)-(23) may be used in order to determine the three unknown 
quantities |n;x 1 ,x 2 ), E n and N n (yi, y 2 ; pi, p 2 ) recursively using quantum-mechanical per- 
turbation theory around the static solution. For this purpose a convenient way to rewrite 
Eqs. (21)-(23) is (E n (x;p) = £ n ( Xl , x 2 ; Pl , p 2 ) and E^(x) = E n ( Xl , x 2 )): 2 

iVn(yi, y 2 ; pi, P2 )5 (3) (yi - xi)5 (3) (y 2 - x 2 ) = 5 (3) (yi - xi)6®(y 2 - x 2 ) 

/V fr 3 (n; yi,y 2 |m; Zl ,z 2 ) ( ) 



E%>(z)-EZ"(x) 

{ Jd 3 x' 2 (°)<m; z 1? z 2 |n; x' x , x' 2 )(E n (x';p') - E^(x'))S^(^ - Xl )5 (3) (x 2 - x 2 ) 

-(°\m-,z 1 ,z 2 \H I \n-,x 1 ,XL 2 )}, (24) 



X 



2 A slightly different set of equations can be found in Ref. [15]. 
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|n;x 1 ,x 2 ) = JcP Zl y"rf 3 z 2 |n;z 1 ,z 2 )(°)Ar T 1 / 2 (z 1 ,z 2 ;p 1 ,p 2 )5( 3 )(z 1 -x 1 )5( 3 )(z 2 -x 2 ) 

4- V fd 3 ? !d 3 7 |m;zi,z 2 )(°) 
m+ J J E y m '{z) -E K n '(x) 

x { jd 3 ^ Jd 3 x' 2 (°><m; z 1; z 2 |n; x' 1; x 2 )(£ n (*'; p') - £M(^))<j(3) (x / - Xl )5 (3) (x 2 - x 2 ) 

-^(mjz^z^lmx^x,)}, (25) 

jd 3 x' 1 jd 3 x' 2 N 1 J 2 (y 1 ,y 2 ;p 1 ,p 2 )6^\y 1 - ^ 3 \y 2 - ^ 2 )E n {x'-p')8^{^ - Xl )5 (3) (x 2 - x 2 ) 

= ^ ) (y)^ /2 (yi,y 2 ;Pi,P2)^(yi -x 1 )5 (3) (y 2 -x 2 ) + (0) (n;yi,y 2 |^|n;x 1 ,x 2 ). (26) 

By means of the above equations it is formally possible to obtain, within the framework of 
a 1/m expansion, the "energies" and the "states" of any excitation of the NRQCD Hamil- 
tonian. 

Up to O(Hf), the energy of a generic state labelled n is given by: 

E n (y;p)5^( yi - ^ 3 \y 2 - x 2 ) = E^(y)5^( yi - ^ 3 \y 2 - x 2 ) 
+ (0) (n;yi,y2|^/|n;x 1 ,x 2 )(°) 

- \ Y. Jd 3 z 1 d 3 z 2 (°\n-y 1 ,y 2 \H I \k;z 1 ,z 2 )W(°\^ 



k=£n ' 



X ^ + 



E^(z)-E^y) E%\z)-E%\x) 

l -Yl jd' Zl d 3 z 2 |d 3 6^2 (0) (n;yi,y 2 |^|k; Zl ,z 2 )( ) 

x(°)(k;z 1 ,z 2 |i/ / |n;| 1 ,| 2 )( ) ( )(n; ^, i 2 \H t \ m x 1; x 2 ><°) 



Ef\z) - E&%) Ef\z) - E&%) 

l -Y^ Jd 3 Zl d 3 z 2 y"d 3 6^2 (0) (n;yi,y2|^/|n;^i,C 2 

x (0) (n; €i, ^l^/lk; zx, z 2 )<°> <°> (k; z 1; z 2 |# 7 |n; x 1; x 2 )(°) 



2 



, 2 ; (0) 



^(,)-^( y )4°)(,)-^(0 
+ 1 £ |^i^2/d 3 6rf 3 6 (0) (n;yi,y2|^|k';^,£ 2 )W 

k,k'^n 

x (0) <k'; £ 2 |tf 7 |k; Zl , z 2 )(°) (°)(k; Zl , za^n; Xl , x 2 >W 



* ' *f (*) - i#>(y) £#>(0 - M 0) (l/) + Ef\z) - E$\x) E$\0 - E&(x) 
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+ 0(Hf). (27) 

The expansion up to O(Hj) was considered in [14] in order to obtain the 1/m potential. The 
O(Hj) term was obtained in [15]. The O(Hf) expression is new. A detailed derivation of 
Eq. (27) will be given in Sec. HID. 



C. The NRQCD states at 0(l/m 2 ) 

The states can also be formally expanded in 1/m: 

|n;x!,x 2 ) = |n;xi,x 2 ) (0) + — |n; x 1; x 2 )« + ^-|n; Xl , x 2 ) (2) + • • • . (28) 



m ' m 2 



It is convenient to write the above states in terms of some new states |n;xi,x 2 ), defined 
recursively as (see also Ref. [15]): 

|n;x 1; x 2 > = Imx^x^ ) + 1 ^ j d^d^m; x^, x' 2 )( '(°»(m; x^, x 2 | 

E n (x) - Hy > m _, n J 

x Xi, x 2 ) — J d 3 x' 1 d 3 x' 2 \n;, x' 1; x 2 )(°)(n; x' l5 x 2 |if/|n; x 1; x 2 )| 

= |n;x 1 ,x 2 )(°) + -|n;x 1 ,x 2 )( 1 ) + ^-|n;x 1 ,x 2 )( 2 ) + -- - . (29) 
m m 2 

As a consequence of Eq. (29), it holds that 

(0) (n;xi,x 2 |n; yi ,y 2 ) = 5 {3) (xi - yi)5 (3) (x 2 - y 2 ) (30) 

or equivalently (this equation will become crucial in later sections to simplify some calcula- 
tions) 

(0) (n;x 1 ,x 2 |n;y 1 ,y 2 )« = V* ^ . (31) 
At (9 (1/m), we obtain 

|n; Xl ,x 2 )« = |l;W> = -£ fa d^fe «„ s 2 ><°> ( ° >fe ^? Hm] %*\ ^ ■ (32) 

k+J El'(z) - E n '{x) 

At C(l/m 2 ), we obtain 

|n;xi,x 2 ) (2) = |n;x 1)X2 ) (2) + |n; x 1; x 2 )$ rm , (33) 

where 

\(0) 



X J d %d% 



. , v (0) (k; z 1; z 2 |gW|n; ^ ggff (°>(n; ^ x 1; x 2 )(°> 

(^W-M 0) (x))(^W-^ 0, (0) 



( ^(,)-^(x)X(o-^(,)) J' (34) 
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and the second term, due to the normalization of the state, reads (note that N = 1 + 
/m 2 + ... is Hermitian): 

|n; Xl , x 2 )l 2 ) rm . = -lJ d 3 x[d 3 x' 2 \0; x^TvfVi, x 2 ; pi, p' 2 )S^fa - Xl )5 (3) (x' 2 - x 2 ) 

= - Jd 3 z 1 d 3 z 2 \n;z 1: z 2 )^ 

v V- Lh *c (0) (n; z 1} z 2 |gW|k; ^ g 2 )(°> W (k; ^ ^gW |n; x 1; x 2 )W 
hJ (Ef(0-E^ X ))(Ef(0-E^z)) • (35) 

By using Eq. (15) and the identities obtained in Refs. [14,15], explicit expressions for 
the above Eqs. (32) and (33) can be worked out. In particular, at order 1/mwe obtain (the 
spin-independent part was first obtained in [14]): 

|n;xi,x 2 ) (1) = 

feW (^)_4o ))2 ( ^)_ £; (0) )2( ^)_ £ (0) ) 
+ 2(Vl ^ ) '(M 0) -4° ) )3 +Vl -(^-^) 2 

, c F (^(fc^BxIn) 

+ [^Ei -> , ^Bx -> , o-! -> <r 2 , Vi -> V 2 , D x -> D c2 ], (36) 

where \n)^ stands for a shorthand notation of |n; Xi, x 2 )(°\ the state that encodes the 
gluonic content of the state |n; x 1; x 2 )(°) and is normalized as ^{n\m)^ = 5 nm (for a precise 
definition, see Eq. (53) and the following discussion). We will use expression (36) in the 
subsequent sections. 

D. Im.Eo with relative accuracy 0{l/m 2 ): structure of the calculation 

In this paper, we are interested in computing ImE^ (actually lmE ) with relative ac- 
curacy 0{l/m 2 ). We will now explain in detail how the different terms of Eq. (27) appear 
within the quantum-mechanical calculation. 

Equations. (24)-(26), as well as the analogous equations in Ref. [15], implicitly assume 
that the Hamiltonian is Hermitian. This is not true at arbitrary orders and the iteration 
of imaginary-dependent terms may lead to problems. Nevertheless, at the relative 0(l/m 2 ) 
accuracy we are aiming at in this paper for the imaginary terms and for the n = state, 
such effects are zero. Therefore, effectively, we only have to compute the expectation value 
of the imaginary part of the NRQCD Hamiltonian in terms of the 0(l/m 2 ) eigenstates of 
the Hermitian part of the NRQCD Hamiltonian. 3 The reason is that the only imaginary 




3 However, a systematic method to work with unstable particles should be worked out if a higher 
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contribution to the states up to 0(1/ m 2 ) comes from the first line of Eq. (34) and this term 
is zero for n = because of the subsequent Eq. (69). 

The imaginary terms in the NRQCD Lagrangian only appear in the matching coefficients 
of the 4-fermion operators, i.e. in £4-/. Therefore, the imaginary part of the NRQCD Hamil- 
tonian has the structure of Eq. (16). Profiting from this structure of the imaginary terms 
and since the iteration of the leading imaginary terms gives zero, lmE can be computed 
from 

Im£ 6 (3) (xi -xi)<J (3) (x 2 -x' 2 ) = (0;xi,x 2 |lm if |0;xi,34). (37) 

Expanding in 1/m the states and ImH, we can identify the different terms of ImEo in the 
1/m expansion: 

Im E = -^Im E [ i ] + -^Im E (3) + -^Im E^ + • • • . (38) 
m z rrr m 4 

They read 

ImEfc^X! -xi)<^(x 2 -x' 2 ) = (0) (0;x 1 ,x 2 |Imif( 2 )|0;x / 1 ,x 2 )W, (39) 



Im Ef 5 {3) (xi - x[)5 {3) (x 2 - x' 2 ) (40) 
= {1) (0;xi,x 2 |lm | 0; xi , x 2 ) (°) + (°) (0; Xl , x 2 1 Im | 0; xi , x 2 ) « , 

Im^ 4) 6® ( Xl - xi)5( 3 )(x 2 - x 2 ) (41) 
= (0) (0; Xl , x 2 |Im |0; xi, x' 2 ) <°> + « (0; Xl , x 2 |Im #( 2 ) |0; xi, x' 2 )« 
+( 2 )(0; Xl , x 2 |Im (2) |0; xi, x 2 )W + ( )(0; x 1; x 2 |Imif ( 2 )|0; x' 1? x' 2 >( 2 ) 
+(°)(0; Xl , x 2 |lmif( 2 )|0; xi, x 2 )£) rm . + £L(0; x 1? x 2 |lmif( 2 )|0; xi, x 2 >^. 

After an explicit calculation, we have 

lm^ 3) =0, (42) 

since 

W(0; Xl , x 2 |lmij( 2 )|0; xi, x 2 )W = ^(0; x 1; x 2 |lm#( 2 )|0; xi, x' 2 )« = . (43) 
Moreover, we have 

( 2 )(0; Xl , x 2 |lmij( 2 )|0; xi, x 2 )^ = (°>(0; x x , x 2 |lmi/( 2 )|0; xi, x 2 )( 2 ) = . (44) 

These results follow from Eq. (31), supplemented by the following argument. The colour 
structure of ImH^ is such that, at the gluonic level, the following matrix elements are 
produced within the total expression: 



precision is warranted. 
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(%|li<g>l 2 |0)( > = ( ><ra|0)( >=<S, 



(45) 



(by definition) and 

W(n\T?®T} a \0)W. (46) 

In order to deal with this second expression, we note that the lowest excitation, in the limit 
xi — > x 2 , has no gluonic content and behaves like |0;xi,x 2 )^ = l c / \/W c \v&c) , so that: 

<°><n|7? ® r 2 ta |0)(°)^ 3 )( Xl - x 2 ) = C7 /< 5 n0 (5W(x 1 - x 2 ), (47) 

where C/ = (iV 2 — l)/(2iV c ). The above expressions may appear problematic since they 
involve the behaviour of the state in the limit xi — > x 2 . and some regularization could be 
required in this case. However, we actually only need a weaker condition to ensure that Eq. 
(44) is zero. What we have is an expression like 

(0) ( 1 Oi | n) {0) ( • • • ) (0) ( A: | ® T 2 ta 1 ) (0) <5 (3) (x! - x 2 ) , (48) 

n^O k^O 

where 0\ is some unspecified operator. Following Ref. [14], this expression is the spectral 
decomposition of the Wilson loop (for the definition of a Wilson loop with a number n of 
operator insertions, see Ref. [15]): 

j dh ... dt n «Oi(ti)(. • -)T? <g> T\ a (t n ))) c <f( 3 >( Xl - x 2 ), (49) 

where ((O)) stays for the insertion of the operator O on a static Wilson loop of spatial 
extension xi — x 2 . In the presence of more operators, the symbol ((• • -)) c indicates the 
connected part (see in particular the erratum of Ref. [15]). One can see that the operator 
(49) is zero in the limit xi — > x 2 . In order to obtain this result, it is very important that 
the delta acts directly on the states 4 . In this situation, and in the limit x x — > x 2 , one can 
see that the disconnected piece of the Wilson loop cancels with the connected piece, proving 
Eq. (44). 

For the other terms, we have 

( 1 )(0;x 1 ,x 2 |lm^ 2 )|0;x' 1 ,x 2 )( 1 ) = 

x (0) (k^ 1 ,| 2 |ImijfJ / |k;z 1 ,z 2 )( )( )(k;z 1 ,z 2 |^ 1 )|0;x 1 ,x 2 )( ) 



4 We may have situations where the Wilson-loop operator has the structure 

J dti...dt n (vi((Oi(ti)(--0TT®r 2 ta (t n ))) c )^ 3 )(x 1 -x 2 ). (50) 

In this case the argument does not apply since the delta does not act directly on the Wilson loop. 
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1 1 1 1 



<°> (0; Xl , x 2 |Im |0; x' x , x 2 )( 2 ) m . + gL (0; Xl , x 2 |Im |0; x' x , x' 2 > <°) 



(0) 



x( )(k;z 1 ,z 2 |//W|0;^ 1 ,O (0)(0 ^0;£ 1 ^ 2 |Im//i 2 _ ) / |0;x 1 ,x 2 \( ) 



x- 



Jd 3 zidh 2 |^irf 3 6 (0) (0;yi,y2|im//| 2 J / |0;^ 1 ,O (0) 

x (0) (0; £ 2 |# (1) |k; Zl , z 2 )(°) W (k; Zl , z 2 |//« |0; x 1; x 2 )(°) 
1 1 

^r^)-4 0) (y)4 0) W-^ u; (0 



X „(0)/ X p(0)/„ X p(0)/ X 7^(0),- ■ ^ 



Indeed, the last two equations hold as well for an arbitrary n and not only for the state n — 0, 
for which we have explicitly displayed them. It can be easily checked that the imaginary 
part of Eq. (27) for n = coincides with the above expression (38) supplemented by Eqs. 
(39)-(41), (43), (44), (51) and (52). 

E. ImE'o with relative accuracy 0(l/m 2 ): explicit expressions in terms of gluonic fields 

The expressions obtained in the previous section can be rearranged in terms of the pure 
gluonic content (see Refs. [14,15]). In order to achieve this we have to make the quark field 
content of the states explicit and use the Wick theorem. There is some freedom in choosing 
the specific realization of the quark fields under spin transformations. In [14], the following 
state was chosen 

|n;xi,x 2 >(°> =^ t (x 1 )x(x 2 )|n;x 1 ,x 2 )W Vxi,x 2 . (53) 

In the basis of 4-fermion operators that we are using in this paper (see Appendix A) and 
in the above basis, the quantum-mechanical operators that naturally appear are l s <g> l s and 
a % ® er 7 , where \ s (<r l ) is the identity (sigma-matrix) in spin space acting either on the 
final- or the initial- spin quark-antiquark state. Analogous definitions can be made for the 
operators acting on the colour subspace. 
Another possibility is the state 

|n; Xl ,x 2 )(°) = V t (x 1 )xJ(x 2 )|n;x 1 ,x 2 )(°) V Xl ,x 2 , (54) 

which has been used in Ref. [15]. The quantum-mechanical operators, which naturally appear 
in this way, are li j2 , cr\ 2 , and they represent the operators acting either on the particle 1 
or 2 (in this case we have always a particle interpretation). Analogous definitions can be 
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made for the operators acting on the colour subspace. This representation appears to be 
more convenient for the calculations of the quantum-mechanical matching. In principle, one 
could also write the local 4-fermion operators in a basis convenient for these states by using 
Fierz transformations [20]. 

In both cases, we assume the state to be properly normalized in the spin sector. Depend- 
ing on the calculation, one definition turns out to be more useful than the other. In any 
case, at the end, we are interested to write the quantum-mechanical Hamiltonian relevant 
to the Schrodinger equation. A way of avoiding ambiguities is to write everything in terms 
of a definite set of spin operators. We will adopt the operators S* and 1 acting on a generic 
1/2 ® 1/2 spin space and defined as 

S l = ^ ® h + li ® Y' 1 = 11 ® h (55) 
It is possible to transform them in the operators l s <g> l s and cr % <g> cr J by using the identities: 
XcS 1 S»"xJ = XV ® <r* X , Xc (2 1 - S 2 ) X J = ^ g, lsX . ( 56 ) 

Let us now compute the different matrix elements that appear in Eq. (41). The contri- 
bution due to the dimension-8 4-fermion operators reads 

(0 H0-,^2\lmH^ f \0-, yi ,y 2 )^ = (c A lm /i( M+1 Pj)7£V\5< 3 >(r) V* 

+ ^Img^Sj^jfrv* + ^i,^(r)} + ^Im f^Pj)r^\v)^ 



x5( 3 )(x 1 -y 1 )5 (3) (x 2 -y 2 ), (57) 

where Ca = N c , V = V r , r = x x — x 2 and (7^ will be used in Sec. V) 

73? = * y (21-S 2 ), (58) 

7^ = ^, (59) 

= 2 €kieek i e ' ^ ^ ' (^0) 

^ = ^ fc s^ + ^s fc _ ^ fc s^ + ^s fc _ (gi) 

nS, = * y (21-S 2 ), (62) 

^ = ^S 2 , (63) 

fi? 1 (^ 1 , 3 £> 1 ) = S i S'"-^S 2 , (64) 

= |flgtf. (65) 



Equations. (58)-(61) and (65) provide the explicit expressions of the operators Ts and T^j 
first used in Ref. [16]. The non-perturbative constant S\ (as well as all the other constants 
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appearing in this section) will be defined in Sec. Ill F. If we 
consider the electromagnetic contribution due to H^f, we obtain (in this case there are no 
octet operators) 



(0) 



(0; xi , x 2 1 Im H^ M) |Q; y i , y 2 



C A Im/ EM ( 2S+1 Pj)T^V^( 3 )(r)V^ 



(66) 



S^jv^ + f ^i,5 (3) (r)} ]5 (3) (xi - yi)5 (3) (x 2 - y 2 ). 



In order to calculate the contribution due to the 1/m correction to the state, we need to 
know (a 1 is understood where no spin-operator is displayed): 



1 (°)(n|[Di.,0Ei]|m)(°) 



E {0) - E {0) 



(o). 



(o) . (o) 



(j\gEi\m) 



(o) 



(°)/r7lnF,J™\(°) (0)/Ti.l/iF.,lm\(°) 



ynin J-Jm ) l- J ri -E/m 

■^cn • ^(nbBilm)^) ^(x! - yi )5 (3) (x 2 - y 2 



^ , gB 1 -> o"! -> <x 2 , Vi -> V 2 , D x 



D 



c2 



Vn 7^ m, (67) 



W(n; Xl ,x 2 |Irn//| 2 _ ) / |m;y 1 ,y 2 )(°) = 
- ( 2^Im/ 1 ( 1 5 )-^Im/ 8 ( 1 ^ 

+ (im^S,) -Imf.CSo) + ^ (lm/ 8 (%) -Im/sf^))^ S 

+ plm/sC^o) + (lm/ 8 ( 3 ^) - Im/g^o)) S 2 ] 5^(r)T F 5 n 
x^Cxx-yO^Cxa-ya), 



5 (3) (r)(n|l c <g> l c |m) 



(68) 



where Fj = F( Xj ), F being a generic gluonic operator. In particular from the last equation 
it follows that 



(0), 



n^Xallmifj^Ojy^ya)* ) = Vn ^ 0. 



(69) 



It is this equation that guarantees that, for the n = quarkonium state, no imaginary 
contribution is carried by the state (see the discussion at the beginning of Sec. HID). 
Finally, from the above equations it follows that the contributions due to the 1/m correction 
to the state read: 

^(0;^ 1 ,^ 2 \lmH^\0;y 1 ,y 2 )^ ^ 
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S 3 V5^(r)V [4Im/ 8 fS ) - 2S 2 (im/sO^o) - Im/s^i))] 
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+ 2T F c 2 F B 1 5 (3) (r) 



lmf 8 ( s S 1 ) + ^S 2 (lmf 8 CS ) - 3Im/ 8 ( 3 1 S 1 )) 



+ ^3 (2)5(3) (r) l41mf 8 CS ) - 2S 2 (lm/ 8 (^o) - Im/ 8 ( 3 ^))] 

- ^(4 2) " 4 2 ' C V (3) (r) [4ImA(^ ) - 2S 2 (im^So) - Im/^ 3 ^))] 
x^( Xl - yi )^(x 2 -y 2 ). 



(70) 



For the electromagnetic contribution we have the intermediate vacuum, which does not allow 
an intermediate emission of gluons. This means that 



( 1 )(0;x 1 ,x 2 |Imi7^|0;y 1 ,y 2 )« = 0. 
The contributions due to the normalization of the state read 
<°> (0; Xl , x 2 |Im tf( 2 ) |0; yi , y 2 ) r ( 1 2 ) rm . + r (2) ra , (0; x 1; x 2 |Im tf (2) |0; y 1; y 2 ) (0) = 



(71) 



(-^ 3 {V 2 ,^(r)} 



Im/xOSo) + y (Wi( 3 Si) - Im/iO^o)) 



-2C A 4^5( 3 )(r) 



Im/xC^o) + y (Im/x^x) - 31™/^%)) 



) ^)(r) [4ImA(^o) - 2S 2 (im/^ 1 ^) - WiW)] ) 
x^( Xl - yi )^(x 2 -y 2 ). (72) 

Exactly the same contribution is obtained from the electromagnetic terms if we change the 
subscript 1 in the matching coefficients by EM. 



F. Gluonic correlators 

The non-perturbative constants £ n , B n , Sf \ S^ 2 '^ and ( e^ 2 ' norm -) ) w hich appeared in the 
previous section, are pure gluonic quantities, since the fermionic fields have been integrated 
out. Within the quantum-mechanical matching, they are first obtained in terms of gluonic 
states. For instance, we obtain the expressions 



f 5 Ji 

n 3 
B ^ 



,(2,c) 



+1 , ^ (%E^K%g|0) 
^ (E (0) - £ (0 V+i ' 



i) n+l n\ 



(0|<7B^)(%B^|0> 



\ (0|^Ei|r> • (r|^Ei|n>(n|^Ei|s) • (s|^Ei|0> 



,(0)' 



n+l 



4 ,e„ I (4 0) - ^° W - ^ 0) ) 4 (4 0) - ^ (0) ) 

(0|^|r) • (r\gEl\n)(n\gB^\s) ■ (s\gE^\0) 



^(O)wp(O) 



(oh 



,(0) 



(0)> 



(73) 
(74) 



17 



(OjgEijr) ■ (rlgE^jnlgE^s) ■ (4gEf]0) 
(4 0) -^ 0) )(4 0) -^ 0) ) 4 (4 0) -^ (0) ) 

(4 0) -^ 0) )(4 0) -^ 0) ) 4 (4 0) -^ 0) ) 



(75) 



For the first two equations, there is no need to specify whether the gluonic fields are inserted 
on the particle or on the antiparticle line since they give the same contribution. We do not 
give here the complete list of expressions at the quantum-mechanical level, since this section 
does this in terms of Wilson loop operators. The former may be derived straightforwardly 
from the latter by spectral decomposition. 

Using the techniques of Refs. [14,15], it is possible to express £ n , B n , 8%'^ and 
£p,norm.) ^ erms Q f ^he more familiar gluonic field correlators. We obtain (traces as 
well as suitable Schwinger lines connecting the gluon fields are understood if not explic- 
itly displayed) 5 

1 f°° 

Sn= N~cJ dttn (9V(t)-gV(0)), (76) 
B n = ±-J™dtt n (gB(t)-gB(0)), (77) 
1 f°° Z"' 1 f t2 

£s W Jo dtl I dh J dt ^^-^) 3 ({9Hh);9E(t 2 )}{gE(t 3 );gE(0)}) c , (78) 



4 2 ' norm - } = ~<l d tl l dt 2 I dt 3 



4^ (y 2 

x ( ((* 2 - t 3 f + (t, - t 3 f) {{gEihh gE(t 2 )} {gE(t 3 ); gE(0)}) ( 
+(*! - t 2 ) 3 ({^E i (t 1 ),^'(t 2 )} {gE%),gW(0)}) c 
+4(ti - hfigEXtJgWitz) gW {t 3 )gE\Q)) ^ 

dh J dt 2 {t l -t 2 f( K (gE\t l )\iT>\gW}(t 2 )gW(<d)) 

+(gE i (t 1 )gW(t 2 )\iB\gW}(0)) + (gE%)\iW , gW](t 2 )gE\0)) ^ 

poo 111 

+ J dt 1 tl(gE%t 1 )[ i B\[ l B\gW]](0))[ + -S S i + -S^, (79) 



3 Note that the quantity £ used in Ref. [16] corresponds here to N c £ 3 . 
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£ ? ] = T^ dt ^\ ^3(t2-t3) 3 {({^ 1 )-^E(t 2 )}{^E(t 3 )-^E(0)}) c 

- A^r^E^) • ^E(i 2 )) Tr(gE(t 3 ) • <?E(0))) C }, (80) 

where 

{{gEfo); gE(t 2 )} {gE(t 3 ) ■ gE(0)}) c = ({gE^) ■ gE(t 2 )} {gE(t 3 ) ■ gE(0)}) 

~ (jE(ii) • gE(t 2 ))(gE(t 3 ) ■ gE(0)), (81) 

and similarly for the other structures with four chromoelectric fields that appear in Eqs. 
(79) and (80). 

For further use, we also define: 

+e 3 [ 2 ' noim -\ (82) 

^2,EM) _ g(2,c) + ^2,norm.) 



IV. PNRQCD 
A. Matching to pNRQCD 

Expressions (27) and alike are no more than formal expansions in Hj, i.e. in 1/m, until 
some dynamical assumption is made. We will assume a mass gap of order A QCD ^> mv 2 
between the lowest-lying excitation and the higher ones. Under this assumption all the 
excitations (n ^ 0) decouple from the ground state (n = 0), which is identified as the 
only degree of freedom of pNRQCD. It corresponds to the singlet state S in the pNRQCD 
Lagrangian (1). Moreover, the above expansion acquires a dynamical meaning, becoming an 
expansion in Aqcd/^ and v in the effective field theory. 

The above assumption is the same as was made in Refs. [14,15] in the situation without 
massless fermions. In this work, we are including light fermions. Nevertheless, at least 
in this paper, we will assume that this does not change the structure of the leading order 
solution (this was also assumed in Ref. [16]). In other words, we will assume that the 
size of the typical splittings between the ground state (heavy quarkonium) and the gluonic 
excitations (hybrids) is much larger than the typical splittings produced by the solutions of 
the Schrodinger equation for the heavy quarkonium. This is, indeed, supported by lattice 
simulations where the plots of the static potentials for the heavy quarkonium and hybrids 
show the same pattern after the inclusion of light fermions [21]. Nevertheless, in principle, 
a new problem may arise. Once light fermions have been incorporated into the spectrum, 
new gauge-invariant states appear besides the heavy quarkonium, hybrids and glueballs. On 
the one hand, we have the states with no heavy quark content. Due to chiral symmetry, 
there is a mass gap, of 0(A X ), between the Goldstone bosons, which are massless in the 
chiral limit, and the rest of the spectrum. We will consider that the Goldstone bosons are 
ultrasoft degrees of freedom and that A x ~ Aqcd, so that the rest of the spectrum should 
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be integrated out. Besides these, we also have bound states made of one heavy quark and 
light quarks. In practice, we are considering the Qq-Qq system. The energy of this system 
is, according to the HQET counting rules [22] : 

m Qg + m Qq = 2m + 2 A. (84) 

Therefore, since A ~ Aqcd, we will assume that they also have to be integrated out. Problems 
may appear if we try to study the heavy quarkonium near threshold. In this case there is 
no mass gap between the heavy quarkonium and the creation of a Qq-Qq pair. Thus, if 
we want to study the heavy quarkonium near threshold, we should include these degrees of 
freedom in the spectrum (for a model-dependent approach to this situation see, for instance, 
[23]). We will not do so in this paper. It may happen, however, that the mixing between 
the heavy quarkonium and the Qq-Qq is small. Indeed, such a mixing is suppressed in the 
large A" c counting. 

Summarizing, light fermions contribute within this picture in three ways: 

1) hard light fermions, they are encoded into the matching coefficients of the NRQCD La- 
grangian and obtained from the computation of perturbative Feynman diagrams at the scale 
to; 

2) soft light fermions, a term that denotes, in a generic way, all the fermions that are 
incorporated in the potentials; it is expected that their main effects can be simulated by a 
variation of the value of the parameters in the potentials; 

3) ultrasoft light fermions, these are the ones that will become pions and, since they are 
also ultrasoft degrees of freedom, they should be incorporated in the effective Lagrangian 
together with the heavy quarkonium. However, we will not consider them in the present 
paper, even if we do not expect to find conceptual problems in an eventual incorporation. 

In conclusion, the matching condition to pNRQCD for the real part reads 

ReE = Reh = + V {0) + + —5- + • • • . (85) 

TO TO TO 

At 0(l/m) the matching has been performed in Ref. [14] and at 0(l/m 2 ) in Ref. [15] (for the 
case without light fermions). We refer to those articles for further details about the structure 
of the potentials. For the imaginary piece, we have the analogous matching condition: 

Im E = Im h . (86) 

Using the results of the previous sections, we can now write the first two terms in the 1/to 
expansion of Imh (the P- wave-dependent terms were obtained in Ref. [16]): 

l mh = + — + ••., 87 

TO TO 

where 

l m/l (2) = _^(3) (r) ^ 4Im/l (i 5o) _ 2S 2 (i m/l (iSo) -Imf^Sj)) 

+4Im/ EM ( 1 S ) -2S 2 (Im/E M (%) -Ihi/em^O) 
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Imh^ = C A T£V^ 3 \r)Vi (im h{ 2S+1 Pj) + Im/ EM ( 25+1 P/)) (89) 

+ ^!,^(r)} (lm gi C s+1 Sj) + lm gEM e S+1 Sj)) 
+ ^T 5 ^ 1 ^)(r)Im/ 8 ( 2 ^P J ) 

+ ^ 3 V5( 3 )(r)V^4Im/ 8 ( 1 5 ) - 2 S 2 (im - Im / 8 ( 3 ^)) 

+2 7>4B! <j( 3 )(r)^Im/ 8 ( 3 S' 1 ) + ^S 2 (Im/ 8 ( 1 ( 5 ) - 3lmf 8 ( 3 S 1 )) 
+|4V(r)(4Im/ 8 (%)-2S 2 (lm/ 8 (%) - Im/^)) 
_^^(3) (r) ^ Im/l( i 5o) _2S 2 (imAC^o) ~Imf 1 CS 1 )) 
-C A ls 3 {V 2 ,5^(r)} (lmf 1 ( 1 S )+lmf EM ( 1 S ) 



+y (Im/i^) -Im/iC^o) +Im/ EM ( 3 5 1 ) -Im/EMC^o)) 
-2 C A 4 Si ^ (3) (r) (im /i(^o) + Im /em(^o) 

+ y (im/i^i) - 31mA (^o) + Im/ EM ( 3 S'i) - 3Im/ EM ( 1 -5 )) 
-^4 2 ' EM) ^ (3) (r) (4 Im huCSo) - 2 S 2 (im huCS ) - Im / EM ( 3 Si)) ) . 

The above expressions have been given in 4 dimensions. Therefore, they should be 
generalized to d dimensions if we want to work in an MS-like scheme in order to use the 
same scheme as for the NRQCD matching coefficients computation. This becomes relevant 
when logarithmic ultraviolet divergences appear in the non-perturbative constants. Hence, 
eventual lattice calculations must be converted to MS in this case. Nevertheless, in several 
situations, it is not necessary to work in an MS scheme if we only want to obtain the 
non-perturbative objects from experiment, since the scheme dependence simply goes into 
a redefinition of the non-perturbative constants. Finally, note also that in addition to the 
divergences in the non-perturbative constants, which are due to large momentum transfers 
k, at some point there will also be ultraviolet divergences arising in quantum-mechanical 
perturbation theory, which are due to large relative momenta p. These must also be regulated 
in dimensional regularization and MS-subtracted, along the lines worked out in Ref . [24] . 



B. Power counting in pNRQCD 

With the above results, we are in a position to compute the inclusive decays of heavy 
quarkonium into light particles by using Eq. (2). Before doing so, we have to specify 
some power-counting rules in order to estimate the importance of the different terms of the 
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pNRQCD Hamiltonian. Previous discussions on this subject, some of which we will repeat 
here, can be found in Refs. [14,15]. 

With the results of Sec. IV A and using Eq. (2), the decay width of S- wave quarkonium 
has schematically the following structure: 

r^imcf 6 '^ ^' 2 (1 + ^ + ..) 

3 m 2 \ m 2 J 

+ Im cf 8 f f W0)(V 2 JW0)) + IU°)I 2 A QCD +..)+... ^ (90 ) 
3 \ to 4 m 2 m 2 J 

where c 4 _/ stands for the NRQCD 4-fermion matching coefficients and R ns os is the S-wave 
radial component of the solution of the real piece of the Schrodinger equation: 

(Re h) (f>nji s (r) = E nj i s <f> nj i s (r), (91) 

with the normalization (|s) sp in denotes the normalized spin component): 

0n S o s (r) = Rns0s(r)-7=\s) spin . (92) 

V47T 

Although E n jis coincides with the binding energy of the system at the order we are working 
at, it will no longer be so when iterations of imaginary parts start playing a role. 

From Eq. (90), we can see how the power counting has to be organized. On the one 
hand, we have an explicit expansion in Aqcd/^, independent of the details of the bound 
state. In the most conservative situation (Aqcd ~ mv), it would correspond to having the 
power counting Aqcd/^ ~ v. We can also find derivatives of the wave function divided by 
to. They typically scale like V/to ~ v. On the other hand, the normalization condition 
of the wave function sets the scaling \R n ji s \ 2 ~ (mv) 3 . This means that a formal 0(mv 5 ) 
accuracy (leaving aside possible ct s (m) suppressions due to the NRQCD matching coeffi- 
cients) is achieved with Eq. (90). At the same order of accuracy, the decay width of P-wave 
quarkonium has the structure: 

r~Im < j= 8 f |Vi ^ 1 ; (0)|2 + ---- (93) 

3 771 

In the above discussion, we have only considered the leading order power counting of the 
wave function at the origin ~ (mv) 3 . This accuracy is sufficient for the P-wave function 
of Eq. (93), as well as for the wave functions multiplying Aq CD /to 2 terms or with two V 
in Eq. (90) but not for the leading order term. In this case, one has to take into account 
that the wave function at the origin also has subleading contributions in v. \R"njis(^)\ ^ 

(mv) 3 (l +av + bv 2 H ). Therefore, we have to further specify the solution of Eq. (91), for 

which we have to set the power counting of the potentials in the Schrodinger equation. Since 
we do not know the specific dynamics of the different potentials, the only thing we can do is 
to require consistency of the theory and allow, in principle, the most conservative counting. 
This would correspond to setting the counting by the largest scale that has been integrated 
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out, i.e. the potentials would scale like (mv) d , d being their dimension. 6 For definiteness, 
we will also assume a s (m) ~ v q with q > 1. 

Leading order. Consistency of the theory requires the virial theorem to be fulfilled. In 
other words, the potential at leading order needs to fulfil 



with the power counting 



p 2 /m ~ Vlo ~ ^1 ~ mf 2 . (95) 



It follows that ~ mt) 2 (even if, using the most conservative power counting, we would 
have obtained ~ mv). Moreover, in our power counting we have V^/m ~ mv 2 . 7 
Therefore, in the most conservative situation, we would have 

, ^ 

V LO = V® + . (96) 

m 

The important point here is that, at this order, the potential is spin- independent (E^ ls = 
E^i and R^]i s = -^i?)- Therefore, the leading-order P-wave function reads 

0i°i ) s (r)=i?i° 1 ) (r)(f|j S ), (97) 

where |f) is the normalized eigenstate of the position and \js) stay for J (total angular 
momentum) and S eigenstates such that 

(f|jo) = y; m (f)|o) spin (j = / = i), (f|ji) = ^ m (f), (98) 

where m denotes the third component of the angular momentum and detailed expressions 
of yj m (r) can be found in Ref. [26], Appendix B. 

Next-to-leading order. The C(l/m 2 ) potential scales utmost as V^/m 2 ~ mv 3 . 
Therefore, in the most conservative situation, we would have 



6 Notice that our power counting rules are different from those of [1,25]. Whereas ours are meant 
to apply in the situation Aqcd 3> mv 2 , the power counting rules in Refs. [1,25] rather follow the 
counting Aqcd ~ mv 2 . Indeed, if we take Aqcd ~ mv 2 in our results we obtain a similar power 
counting for the NRQCD matrix elements. 

7 As a consequence, if the potential is non-perturbative, we have no general argument to 
consider V^/m subleading with respect to V"(°). A lattice simulation or some model-dependent 
studies are, therefore, highly desirable to discern the issue. Whereas it is difficult to obtain this 
information from the spectrum structure, the study of the decays may perhaps shed some light on 
this problem. Finally, we note that, in the perturbative situation, W 1 ) has an extra a 2 suppression. 
Further discussions can be found in Ref. [14]. 
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\/( 2 ) 

Vnlo = -5-. (99) 



At this order, spin- dependent contributions start to appear. In particular, the spin- 
dependent potential contributing to the S-wave function at the origin reads 



6V=^p-BeV$> 1 \r), (100) 



where [15] 



o '2 poo 

ReV£' 1 \r) = ^ft / dt{(gB 1 (t)-gB 2 (0))) + 2C A (Ref 1 ( 1 S )-Ref 1 ( 3 S 1 )) 5®(r). 
<J Jo 

(101) 

This potential produces the following correction to the S-wave function: 

R ns0s {0) _ B${0) 1 / 3\ 1 ReV^InO^ U02) 

where P n = I — |n0)(n0| and (r\njls) = 0^ s (r) ((r|n0) = 0^°J(r)). 

If the spin-dependent potential (100) is 0(mv 3 ), it just provides the leading order spin- 
dependent correction to the S'-wave function at the origin and one can use the difference 
between vector and pseudoscalar decays to fix the value of the correction. If the spin- 
dependent potential is (9(mt> 4 ), it provides a correction to the S'-wave function squared at 
the origin, which is of the same order as the 0{v 2 ) corrections to the decay width that we 
have already evaluated. Therefore, in this last situation, Eq. (102) would account for the full 
difference between the vector and pseudoscalar wave functions at the origin at relative order 
0(v 2 ), which is the precision we are aiming at in this work. This last counting seems to be 
supported by the size of the spin-dependent splittings in the bottomonium and charmonium 
spectra. 

For the spin-independent contributions, we will make no assumption at this or higher 
orders, as their effects will be encoded into the wave functions, which will be left uneval- 
uated. Our results allow for the most conservative counting where V^/m ~ mv 2 and 
(spin-independent )/m 2 ~ mv 3 . We note that, in this power counting, potentials with 
imaginary part arise in the pNRQCD Hamiltonian at order ma s (m) 2 v 3 (where the powers 
in a s (m) come from the imaginary part of the 4-fermion matching coefficients in NRQCD). 
Therefore, corrections due to the iteration of imaginary terms, which could affect the valid- 
ity of Eq. (2), are far beyond the accuracy of this paper. In fact, the general factorization 
formula put forward in [1] may not hold beyond a certain order. 

In any case, we do not rule out that a different power counting may also lead to consistent 
equations in the non-perturbative regime for some specific ratios of Aqcd versus m and versus 
p and k. This point deserves further investigation and may lead to a different implementation 
of the matching procedure. We recall that the issue of assessing the power counting in the 
non-perturbative situation has been addressed before by Beneke [27] and by Fleming et 
al. [28]. In both cases, the authors have given some freedom to the possible size of the 
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NRQCD matrix elements by introducing a parameter A that interpolates between the power 
counting in the perturbative limit and other possible power countings according to the value 
of A. In this respect, our formalism may shed more light to clarify this problem, since it 
incorporates the factorization between the soft and the ultrasoft scales, allowing us to write 
the NRQCD matrix elements in terms of the wave function at the origin and of some bound- 
state-independent constants. Another point of concern is whether there are non-perturbative 
effects that are not accounted for in the 1/m matching. 

We conclude this section by giving a useful equality, valid in dimensional regularization: 

^ ) (Q)(v 2 4?(Q)) _ |iff(o)l 2 flff (103) 

m 4 m 2 m ' 

which follows from the fact that we know the behaviour of the potential and the wave function 
(up to a constant) at short distances and that (see Appendix D) 

(n, j, I, s\V^\r = 0) = (n, j, I, slV^ |r = 0) = (in dim. regularization). (104) 

With this we have discussed the relative importance of the different terms that will appear 
in our evaluation of the decay widths. The results can be found in Sec. VI. 

V. THE MATCHING IN THE CASE MV > A QC d > MV 2 

Although it is not clear whether quarkonia states fulfilling mv ^> Aqcd fnv 2 exist in 
nature (mv ~ k ~ p and mv 2 ~ E will always be understood in the present section), this 
situation is worth investigating for several reasons. First of all, the calculation in the general 
case of the Sec. Ill is non-standard and, hence, any independent check of it, even if it is in 
a particular case, is welcome. Secondly, the calculation in this case can be divided into two 
steps. The first step can be carried out by a perturbative calculation in a s , which involves 
far more familiar techniques. The second step, even if it is non-perturbative in a s , admits a 
diagrammatic representation, which makes the calculation somewhat more intuitive. Third, 
the more detailed information on the potential allows us to make important tests on how 
the terms in the potential can be consistently reshuffled by means of unitary transformations 
[14], as is illustrated in the example provided in Appendix E. 

A. pNRQCD' 

As mentioned in the Introduction, we shall call pNRQCD' the EFT for energies below 
mv. Since mv ^> Aqcd, the integration of the energy scale mv, namely the matching between 
NRQCD and pNRQCD' can be carried out perturbatively in a s . This is done following Refs. 
[4,6]. A tree-level matching is sufficient, but higher orders in the multipole expansion will 
be needed. We only display below the terms eventually required in the calculation: 

£ pN RQCD' = Tr { St ( id ° ~ h s)S + O f (lD - h Q ) O} 

+Tr I Oh ■ gE S + H.c. + ^ + \ 
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+ ilr {O f rV gB^ O - O+OrV gWW] 
+ 2l Tl { 0tr<rirfe 9 r > i ~ Dj E k S + H.cj 

+ ^Tr {Ot^ - <r 2 ) ■ gB S + H.c.} - ^G"^ , (105) 

where the traces are over colour space only. S and O are chosen here to transform as a 
1/2 ® 1/2 representation in spin space (hence <J\ — cr 2 = cr\ (g) 1 2 — li (g) cr 2 ); h s and /i read 
as follows (again we only display terms eventually required in the calculation): 

h s = -—- C f ^ - i ^-^(umhCSo) - 2S 2 (im A(%) - W^)) 
m r 2 m 2 \ v 

+4Im/ EM ( 1 S , ) -2S 2 (Wem(%) -Wem^i)) 
^^{v'VV^r)} (im^r^^ + Im^Mr^^)) , (106) 

^ = -- + ?^-^)- 
m \ 2 J y r 

.T F 5( 3 )(r 



2 m 2 



4Im/ 8 ( 1 S ) -2 S 2 (lm/ 8 (%) -Im/ 8 ( 3 ^)) 



+i7>7&V— ^V j lmf 8 ( 2S+1 Pj) . (107) 
The Feynman rules associated to this Lagrangian are displayed in Fig. 1. 



B. Matching pNRQCD to pNRQCD' 

The matching of pNRQCD' to pNRQCD can no longer be done perturbatively in a s , but 
it can indeed be done perturbatively in the following ratios of scales: Aqcd/^w (multipole 
expansion), Aq CD /m and mf 2 /Aq CD . The diagrams contributing to the calculation are 
displayed in Figs. 2 to 9. 

We have focused on contributions to S-wave states involving imaginary parts. Since 
the imaginary parts, which are inherited from NRQCD, sit on local (5®(r), V5 ( - 3 - ) (r)V, 
etc.) terms in the pNRQCD' Lagrangian, they tend to cancel when multiplied by the r's 
arising from the multipole expansion. Hence, for an imaginary part to contribute, it is 
necessary to have a sufficient number of derivatives (usually arising from the 171v 2 /Aqcd 
expansion) as to kill all the r's. Since derivatives are always accompanied by powers of 
1/m, it implies that at a given order of 1/m, only a finite number of terms in the multipole 
expansion contribute. In our case a fourth order in the multipole expansion is sufficient. The 
natural way to organize the calculation in our case would be to assign a size mv p to Aqcd, 
1 < p < 2 and v q to a s , 1 < q < 2, and to carry out the calculation at the desired order 
in v. However, our main goal here is not the phenomeno logical relevance of the situation 
mv ^> Aqcd rnv 2 , but providing an independent calculation to support the results of Sec. 
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IV A. Hence, irrespectively of what p and q may be, we will only be interested in fishing up 
the imaginary pieces that contribute to S-wave states up to order 1/m 4 . 

The two diagrams in Fig. 2 correspond to the leading contribution in the Aqcd/^w 
and Aqcd/"^ expansion, respectively. Figure 3 displays the evaluation of each of them in 
the mv 2 /Aqcd expansion. The diagrams in Fig. 4 correspond to the next-to-leading order 
contributions in the KqcT>l mv expansion, and Figs. 5-9 display their evaluation in the 
mf 2 /AqcD expansion. It is then clear that the basic skeleton of the calculation consists of 
the x = (Aq CD /mt>) 2 and y = (Aq CD /m) 2 expansions, which suggests writing the pNRQCD 
Hamiltonian as: 

h = h s + h x + h 2x + h y + .... (108) 
The interpolating fields of pNRQCD' and pNRQCD will be related by: 

•$1 pNRQCD' = Z 5 5'|pNRQCD = (1 + Z x + Z 2x + Z y + ... ) 2 5* | pNRQCD 

= (l + \(z x + Z 2x + Z y - l -Z 2 )j S\ pNRQCD . (109) 

The matching calculation reads: 

J dt e~ iEt J d 3 R (vac|T{ 1 S(R, x, t)S(0, x', 0)} |vac) | p nrqcd' 

= ! dte~ lEt /"rf 3 RZ5( vac |T{ 1 S(R,x,t) 1 S(0,x',0)}|vac)| pNRQ cD^5 t . (110) 

J — oo J 

The right-hand side of the matching calculation has the following structure (up to a global 
% factor, which is dropped): 

1 1 1 1 / „ „ Zl\ 1 



11/ Z 2 

(h X + h 2x + hy) + ~ [ Z X + Z 2x + Zy~^ 



t 



E-h s E — h s y E-h s 2 V y 4 J E - h s 

+ _J_I (z. + z » + z,-2)\ ( z -\ — !— M 

E-h s 2 \ x 2x y A J \2 ) E -h s \2 ) 

+ E^ s hx E^ s K E^ s + {i)E^ s hx E^ ■ (1U) 

Hence, once we have made sure that, up to contact terms, the left-hand side of Eq. (110) has 
exactly this structure, we can easily identify the contributions to the pNRQCD Hamiltonian 
from the second term of the expression (111). 



C. Calculation 

Let us then proceed to the calculation of the left-hand side of Eq. (110) (in order to 
match Eq. (Ill) a global i factor will also be dropped). 
Diagram (a) of Fig. 2 gives: 

— _jf A(lr . 9 E (t )e-'^ lr . g E(0)> — . (112) 
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The fact that mt> 2 /AQCD is small is implemented by expanding the exponential. This guar- 
antees that we will eventually get usual, energy-independent, potentials. 

The first contributions arise at O (?m> 2 /AQ CD ) from the 0{l/m A ) P-wave (Fig. 3a) and 
S-wave (Fig. 3e) terms in the octet potential of (107): 

w ^ ^^ f. tW ), E( o)>|M, (113) 

where are defined in Eqs. (58)-(61) and 7$ in Eq. (65). 

At 0(m 2 v A /AQ CD ) and higher, it is convenient to write E — h Q — E — h s + (V — V s ). Ill- 
defined expressions arise in the calculation, from products of distributions (both products of 
two delta functions and products of delta functions with non-local potentials, which explode 
as r — > 0). It is most convenient to use dimensional regularization in this case, which 
sets all these terms to zero. This is shown in Appendix D, where the relation with other 
regularizations is also discussed. Having this in mind, it is clear that, at the order we are 
interested in, Im (V D — V s ) 2 = and Im (V — V s )r = 0. Hence, we only have to consider: 

v{E-h s fr-*—. (115) 



E-h" E-h 



s 



If we decide to take one power (E — h s ) to the right and one to the left we have: 

r2 + r|r - h - ] Ehr, + E L ir, lh " T]T + Ehr, [h - I]l1 ' '^shr, • (116) 

which does not produce any imaginary part. However, an equally acceptable expression is: 

r2 4l r -l r ' fc -ll¥^ + ^I^ r l- r l + ¥^^l[ r - A J'' l -l' r >^' < 117 > 

which does produce an imaginary part. This apparent paradox only reflects the fact that 
expression (115) by itself (as well as some of the expressions we will find below) does not 
determine uniquely its contribution to the potentials. This expression always leads to contact 
terms, wave-function normalization and potential, as is apparent in (116) and (117), but 
depending on how we decide to organize the calculation, the terms associated to each of 
these pieces change. For instance, when matched to (111), (116) gives: 

h x = [h a ,T][r,h a ] Z x = r[r,h a ], (118) 

whereas (117) gives: 

K = \{[[r, h s ], h s ],r} Z x = i[r, [r, h s }} . (119) 

This should not be a surprise. It has already been discussed in Ref. [14] that this ambi- 
guity exactly corresponds to the freedom of making unitary transformations in a quantum- 
mechanical Hamiltonian, and does not affect physical observables. This is discussed in detail 
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in Appendix E for the decay widths of the S'-wave states we are concerned with. In order 
to fix the contribution to the potential of any term once forever, we will use the following 
prescription. If we have an expression with singlet propagators 1/(E — h s ) only in the exter- 
nal legs, and an even number of powers of (E — h s ), we will take the one closest to the left 
propagator to the left and the one closest to the right propagator to the right, and repeat 
until no power is left except in contact terms. Accordingly, in the intermediate steps, when 
terms with a single external leg 1/(E — h s ) and several powers of (E — h s ) are produced, one 
should take these powers towards the 1/(E — h s ) leg until no power is left except in contact 
terms. If the number of powers of (E — h s ) is odd, we use the same prescription until a 
single power is left. We then write (E — h s ) — (E — h s )/2 + (E — h s )/2 and take one half 
to the right and one half to the left. Expressions with an internal singlet propagator also 
appear, which require a more careful treatment as will be discussed after (128) below. Note 
that this prescription to organize the calculation needs not coincide with the prescription for 
fixing the wave-function normalization in Sec. IV A. Hence, we only expect to agree with 
the results of that section up to a unitary transformation. Anyway, with this prescription, 
(115) gives rise to the potential obtained in (116) and hence to no imaginary part. 

At C>(mV/A| CD ) only the following two terms in (E - h Q f = (E - h s ) 3 + (E- h s )(V s - 
V )(E — h s ) + ... contribute, giving rise to: 

\ {{E - h s )r 2 + v 2 (E- h s )) + I ([h s , r]r + r[r, h.]) 

1/1 1 

+ o p TT \ h 'i \ h *i r ]] r + r ([ r ' h '\> h - 1 



2 \E-h s L ' JJ LL ' J ' 'E-h s 
+\ (W, h a ], h s \ E 1 _ h + E ]_ h [h a , [h s , r]r] 

+ [h s , r] [r, h a ] E \ h + E ]_ h [h s , r] [r, h s \ 

+ \ {e^s [K [K [r ' K] ^h s + ^h s [K r] ks] ' ks] ^ 
+r (V s -V Q )r+ Y^[hs, r] (V, - V a ) r 

+r (V s - V ) [r, h s ]-±— + -et^-tK r] (V s - V Q ) [r, h 1 



E-h s E-h s " iy E — h s 



It is the first term in the fifth line that renders the contribution depicted in Fig. 3d: 



(120) 



At (9(m 4 t> 8 /AQ CD ) and higher, only imaginary parts beyond 1/m 4 are produced. 

Consider next the diagram Fig. 2b. Since the chromomagnetic moment already provides 
two powers of 1/m, only the linear term in the expansion of the exponential contributes 
(Figs. 3b and 3c). This gives: 



r 



E-hs' 
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(c) ^«pM r , i(gB(t , gB(0)) |M. 

Consider next Fig. 4a. Because of the four r's in the expression, only the following term 
in the expansion (E — h Q ) 3 — (E — h s ) 3 + ... contributes. We obtain: 

/s!U!/ir ' L ' S ' 5 W« + Mii + W /""(ft^^E^py.pD* <7E , (0)]]>- <5( ' , ' lr) 



E-h a 72 m 4 " lJ " 1 M J> 1 " J ^J ,l- ^~^nu E _ hs - 

(122) 

For the symmetric diagram, we have: 

E-h. TS ^nm' 1SS) (6ij6kl + 5ik5jl + 5il5jk) r dU3{[B% [D '' » E *W^ E, (°)>^^ • 

(123) 

In fact both contributions are the same, adding up to (see formula i) above Eq. (15) of Ref. 
[14]): 

i TsImf^Ss) /-~ .^un ^mv<5 (3) ( 



E-K 12 m 4 



^°°rftt 3 ([D,^E(t)][D,.^E(0)])|-^. (124) 



Consider next Fig. 4b. The only contributions come from (E — h ) 3 = (E — h s ) 3 + ... in 
one octet propagator and 1 in the other. We obtain (Fig. 5b): 

1 T Im f f 2 ^ +1 S' ) f 00 f tL r 

~B — r S o 4 — -(SijSki + StkSji + 8 u 5 jk ) / dh / dt 2 (ti-h) 3 + t 3 
E - h s 12 J J L 

x^E^tOpy^E^^E'CO))!^. (125) 

Then consider Fig. 4c. From here we get several contributions. Because of the four r's we 
need a total of three powers of (E — h ). When all the powers come from the octet propagator 
in the middle, we get contributions from (E — h ) 3 = (E — h s ) 3 + (E — h s )(V s — V )(E — h s ) + - ■ ■. 
The ones from the second term read (Fig. 6): 



E - h, 6N n m 4 

x 



(T F Im/ 8 ( 2 *+%) - AUmAr+%)) 



poo pti pti 

\ dt x / dt 2 / dt 3 (t 2 -t 3 ) 3 
Jo Jo Jo 

x{({ ff E(t 1 ),.</E(t 2 )}{ ff E(ts),-»E(0)}> 

-±(Tr(gE( tl ) ■ gE(t 2 )) Tr(gE(t 3 ) ■ </E(0))>}|^£ . (126) 

When a power of (E — h a ) does not come from the octet propagator in the middle, all the 
powers can be substituted by (E — h s ). If we put these contributions together with the first 
term before (126), we obtain (Fig. 7): 
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1 



E-K 12m 4 



dh 



dU 



(h - t 3 f + 4 



X 



({ g E(h);gE(h)}{gE(t 3 );gE(0)}) - — (Tr[gE(h) • gE(t 2 )\ Tr[gE(t 3 ) ■ gE(0)]) 



/•oo rti rt-2 

+ dh dh / dt z (h-hf + 4 
Jo Jo Jo 



x 



+ 



{{gE%),g&(t 2 )}{gE%),g&(0)}) - jf (Tr\gE%)gW{t 2 )]Tr\gE%)gW(0)]) 



({gE\h),gW(h)}{gW(h) } gE\0)}} - — (Tr[gE\h)gW (hWgW (t 3 )gE\0)}) 



x 



<S(3)( r ) 

E-K 



(127) 



Consider next Fig. 4d. Clearly this diagram contains the iteration of lower order po- 
tentials, which must be isolated. This is achieved by adding and subtracting the projection 
operator into the gluonic ground state 1 = (1 — |0)(0|) + 10)(0|. The piece (1 — 10)(0|) contains 
new contributions to the potential only, whereas the piece |0)(0| contains both the iteration 
of lower order potentials and new contributions to the potential. Consider first the piece 
(1 — |0)(0|). It is identical to Fig. 4c by taking V Q — > V s in the expression before (126) and 
changing the chromoelectric field correlators accordingly. We then have (Fig. 8): 



1 



E — he 3iV r m 4 



dh 



dh 



dh 



(h - h? + A 



X 



(TrfsE^) • gE{h)\ Tx[gE{t 3 ) • <?E(0)]> - (gE(h) ■ gE(h))(gE(t 3 ) ■ gE(0)) 



/•oo rli pt2 p 

+ / dh I dt 2 I dt 3 (h-hf + t 

Jo Jo Jo L 

x ([(Tr[gE\h)gW(h)]TrigE\h)gW(0)]) - (gE\h) g& (h)) (gV\h) g& (0)) 

(Tr[gE%h)gW(h)} Tr[gW (t 3 ) gE' (0)}) - (gE%)gW(h)) ( 9 W(h)gE\0)) 



x 



6&(r) 
E-h* 



(128) 



Consider next the contribution from |0)(0|. The vacuum insertion leads to an internal 
singlet propagator. To be specific, we have: 



-i 1 
E~=JT a NlJ 



dt (ir ■ gE(t)e 



-i{h -E)t 



ir-^E(O)) 



1 



E - K 



x 



/•oo 

/ dt' (ir ■ gE(t')e-* h °- E)t 'ir ■ gE(Q)) 
Jo 



1 



E - K 



(129) 



The exponentials of (E — h Q ) will be expanded. In order to be consistent with the calculation 
of the lower order potentials and subtract only their iteration, we must treat the powers of 
(E — h ) at each side of the internal singlet propagator exactly as we did in the calculation 
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of the lower order potentials. Let us illustrate how it works when we have two powers of 
(E — h Q ) on each side. The only contributions occur when (E — h ) ~ (E — h s ). If we write 
the propagator in the middle as 1/(E — h s ) = (1/(E — h s )) (E — h s ) (1/(E — h s )) we can use 
(115) and (116) in order to obtain: 

11 1 

r 2 + r[r, h s }- — + -[h s , r]r + —[h a , r][r, h 



E-h s ' E-hs 1 """' 1 ' ' E - /i/~ a '~ jr, '~ SJ £ - h a 
x (E - h s ) (130) 



x 



( r2 + r[r ' hs] Ehs + ^h s [hs > r]r + ^h s [K r] [r? hs 



si E-h s 

We can easily identify the contributions that match the following terms in (111): 
Z x \ 1 f Z x \ ] 1.1.1 



+ T, i- h x-p, i-K 



2jE-h s \2j ' E-hs"* E-h s x E - h s 

+ "(f) E^hs^E^hs^ ' E^hs kx E^hs{^) ■ 

We also see that, apart from the terms above, there are additional terms in (130) that may 
(and do) eventually lead to new contributions to the potential (none of them with imaginary 
parts for this example). For them we use the same prescription as stated at the beginning of 
the section. The contributions to the imaginary parts come from the following terms in the 
expansion only: (i) an (E — h ) 4 from an octet propagator and a 1 from the other one (Fig. 
9, first diagram), and (ii) an (E — h Q ) 3 from an octet propagator and an (E — h ) from the 
other one (Fig. 9, all of them). They read: 



hs-E' 



d£ (<7E(O-<7E(0)> 

1 J 27 N c m 4 h s -E\J Vy V 1 y V ") 

x^dt't'(gE(t')-gE(0)))^^ 



D. Results 

Combining all the above calculations we obtain the same result as in Sec. IV A, except 
for the terms proportional to lm( +1 Ss)- With the mere replacement: 

g (2,EM) _^ ^(2,EM) ^ (132) 
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where we have defined: 



c(2,t) 



x 



dti 



dU 



dU 



({^(txhsE^H^h^O)}) - w (gE(t 1 )-gE(t 2 ))(gE(t 3 )-gE(0)) 



poo pti rt2 r 

+ / dh dt 2 / dt 3 (h-hf + t 
Jo Jo Jo 



X 



+ 



({gE%),gW(t 2 )}{gW(t 3 ),gW(0)}) - — (gE%)gW(t 2 ))(gE%)gW(0)) 
<{^(t 1 ),^'(t 2 )}{^(^),^(0)}) - -l^E i (i 1 )^^ 2 )><^(£ 3 )^E i (0)> 



N r 



poo pt\ 

-KSijSu + Sirfn + SuSjk) dh / dt 2 + 4 (gE%)\jy,gE k ](t 2 )gE l (p)) 

Jo Jo 1 J 

poo 7 2 ^ 

+ j dtt 3 ([D;gE(t)][D;gE(0)]} + -S 4 S + -S 3 S 1 j , 



(133) 



and 



(2,EM) 
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(Tr[^E(t 1 ) • g E(t 2 )]Tr[gE(t 3 ) • <?E(0)]) - (gE^) ■ gE(t 2 )) (gE(t 3 ) ■ gE(0)) 
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Jo Jo Jo L 



x 



({ffE i (O,^'(*2)}{<7E i (t s ),l7E»"(0)}>-— {gWitJgWfofiigWitdgWiO)) 



-i(SijS M + 5 ik Sji + 5 a 5j k ) / dti 



dt 2 



(h - h) 3 + t\ (gE%) py, <?E fc ] (t 2 )^E'(0)) 



+ 



poo 7 2} 

J rftt 3 ([D-^E(t)][D-^E(0)]) + -^o + -^ij , 



(134) 



the same expressions apply. 

As mentioned before, the difference is due to the different prescription to fix the wave- 
function normalization in Sec. Ill B. In Appendix E we show that there exists an unitary 
transformation such that our results can be taken to the form of Sec. IV A, and hence they 
are equivalent for all purposes. 

In fact, it is somewhat surprising that the two calculations lead to identical results (up 
to a unitary transformation). On general grounds, one could only expect that the result in 
this section be a particular case of the general results of Sec. III. In fact the real parts 
of the potentials in h are indeed particular cases of the potentials in [15]. However, since 
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we did not need their specific form at any stage we have not lost generality in our final 
expressions. More surprising is the fact that the matching coefficients of the terms in the 
multipole expansion in pNRQCD' (105) were only calculated at tree level here, whereas the 
expressions in Sec. Ill correspond to an all- loop result. This indicates that there must be a 
symmetry protecting these terms against higher-loop corrections, which may (or may not) 
be an extension of reparametrization invariance [29] or Poincare invariance itself [30]. 8 

In summary, we have presented in this section an alternative derivation of (141)-(146), 
which does not rely so heavily on the 1/m expansion. The matching from NRQCD to 
pNRQCD', which can be done perturbatively in a s , can indeed be implemented in the 1/m 
expansion, as originally proposed [4] , but it can also be done entirely in the framework of the 
threshold expansion [31,13], where the kinetic term is kept in the denominator for potential 
loop contributions and the on-shell condition is used (the results obtained in either way are 
related by local field redefinitions). The matching between pNRQCD' and pNRQCD is done 
in the Aqcn/fnv, Aq CD /m and mf 2 /Aq CD expansions. The approaches taken in these two 
steps are quite different from the strict 1/m expansion of Sec. Ill, and the coincidence of 
the results strongly supports their correctness. 

VI. RESULTS 

In this section we list our expressions for S-wave decays up to 0(c(a s (m))mv 3 x 
(AQ CD /m 2 , E/m)) and for P-wave decays up to 0(c(a s (m))mv 5 ). The P-wave decay widths 
were first obtained in [16] and are given here for completeness. The S*-wave decay widths 
are new. In order to help the reader and for further convenience, we will start by recall- 
ing, at the same level of accuracy, the expressions of the decay widths as they are known 
from NRQCD. In the following we define the radial part of the vector S'-wave function as 
P ral0 i = R^q = P„°o (1 + 0{v)) and the radial part of the pseudoscalar S-wave function as 
P n00 o = P^o = -^no(l + ^( v ))- The quantity R^i' is the derivative of the leading order 
P-wave function. The symbols V and P stand for the vector and pseudoscalar S-wave heavy 
quarkonium and the symbol x f° r the generic P-wave quarkonium (the states x( n 10) and 
x(nJl) are usually called h((n — 1)P) and x.j(( n ~~ 1)-P)> respectively). 

A. Decay Widths in NRQCD 

Including up to the NRQCD 4-fermion operators of dimension 8, the inclusive decays of 
heavy quarkonia are given by: 

T(V Q (nS) -> LH) = ^(lmf l CS 1 )(V Q (nS)\O l ( 3 S l )\V Q (nS)} 
+Imf 8 ( 3 S 1 ){V Q (nS)\0 8 ( 3 S 1 )\V Q (nS))+fo 



For the leading order term, the non-renormalization was verified at one loop in [11]. 
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m 2 m 2 

+Im/>(Jpj («fl»M + Im /tfPi) (W(W) , ( 135 ) 

///- ' ' ///- / 

T(P Q (nS) - LPT) = A Am/ 1 ( 1 5o)(P Q (n5)|0 1 ( 1 5o)|PQ(^)) 

+Im/ 8 ( 1 5 )(PQ(n5)|0 8 ( 1 S'o)|PQ(n5)) + Im/ 8 ( 3 5 1 )(PQ(n 1 S)|0 8 ( 3 1 S 1 )|PQ(^)) 

+Im^%) (M ^ (136) 

m 2 m 2 J 

m 2 y m 2 
+/8( 2S+1 ^)(XQ(nJ5)|0 8 ( 1 5o)|x Q (nJ5))^). (137) 



At the same order the electromagnetic decays are given by: 

2 

m 



r(V<,(nS) -. e + e") = ^ ( Im/ e .( 3 5 1 )<V r <J (n5)|OEM( 3 S,)|V<,(nS)) 



^^QiM^g^MM)), (138) 



T(P Q (nS) -> 77) = Jj( Im/ 77 (%)<P e (nS)|O E M( 1 S )|f > e(nS)) 



+Im977( , So) (M^M(^M^)>Y (189) 

m 2 / 

r(xQ („Jl) -.77) = 2lmf ^ eP /x Q (nJl)\OMy\xdnn)) for = „_ (14Q) 



B. Decay Widths in pNRQCD 

Up to 0(c(a s (m))mv 3 x (Aq CD /m 2 , P/m)) for S-wave and (9(c(a; s (m))7ra; 5 ) for P-wave, 
the inclusive decays of heavy quarkonia are given in pNRQCD by: 



T(V Q (nS)^LH)- CAlR ^ m2 



7T m? 



Im/ 1 ( 3 5 1 )(l-^^ + ^V + C ^ 1 



m 9 3m 2 3m 2 
( 2 ) (n . io _ r*v2 



Im ^3^ 2(^/2-^)^ Imf , (^72-^)4^ 
-Im/ 8 ( 50 — Im/ 8 ( S ) — 2 



+lmg 1 ( 3 S 1 ) 



^nO _ £l 

m m 2 



(lm/ 8 ( 3 P ) + 3Im/ 8 ( 3 P 1 ) + 5Im/ 8 ( 3 P 2 )) ^^J^ 



(141) 
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T(P Q (nS) - LH) 
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-Im/sC^o) 
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Im/ 8 ( 3 ^) 
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m 2 
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r( XQ (nJS) ^ Ltf) = 



r» |p(°)VnM2 



2T f 
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At the same order the electromagnetic decays are given by: 



31m hC b+1 Pj) + £f Im / 8 ( aa+1 S s )53 



(142) 
(143) 



r(V r (nS') -> e+e- 
+Im# ee ( 3 S 1 ) 

r(P Q (n5)^77) 
+Im5f 77 ( 1 S'o) 

r( XQ (nJl)^7 7 ) = 3 



Ca Kq(0)P 
7r m 2 



lmf ee ( 3 S 1 ) 1- 



E$> 2S 3 2£f EM) 



m 9 



+ 



3m 2 3m 2 



(e {0) 




I m 





(144) 



CU|i2S,(0)| : 



7r 77r 



F (0) 9>7 9>7 (2,EM) 2w 

Im/ T7 (^ ) ( 1 - + + ^ 



m 9 



3m 2 







I m 


m z J 



n I pw 



7r ttt 



(°)VnM2 



-/ 77 ( 3 Pj) for J = 0,2 



(145) 
(146) 



C. NRQCD Matrix elements 

By comparing the decay widths in NRQCD and pNRQCD we obtain the following dic- 
tionary between the matrix elements of NRQCD and the non-perturbative constants of 
pNRQCD, valid up to (once normalized to m) 0(v 3 x (AQ CD /m 2 , E/m)) for the S-wave 
matrix elements and up to 0(v 5 ) for P-wave matrix elements: 

KMMSJWnS)) = C A \&>1 (l - <f + ^ + ^) , (147) 

I rn^l I 2 / OF 9 <r(2,EM) 2 » \ 

(V^OMV^nS)) = C ,M_ ^ _ ^ + ?|_ + ^) , (149) 
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\R P (r\\\ 2 ( OP 9 c(2,EM) 2 » \ 

(P (^)|O EM (^ )|P (^)) = ^ _ ^ + + ^ j , (150) 

= ^ipir(o)i 2 , (i5i) 

(VQ(n5)|n( 3 ^i)|Vb(^)) = (P Q (n5)|P 1 ( 1 5 )|P Q (n5)) 

= (Vq^IPemC^OIVq^)) = (Pq^IPemI 1 ^)^^^)) 



(V Q (nS)\0 8 ( 3 S 1 )\V Q (nS)) = (P Q (nS)\O s CS )\P Q (nS)) 



J n0 



£1) , (152) 



r ( 2(C A /2 - C f )4 2) , 

6A 2tt 3^ 1 ' (153) 



= ^^^l 3^ I' (154) 



(x Q (nJ5)|0 8 ( 1 5o)|x Q (^)) = -f ^ l 2 Ujl g 3 - (156) 

Any other S'-wave dimension-6 matrix element is at NNLO and any other S-wave 
dimension-8 matrix element is at LO. 

Equation (152) is worth emphasizing. It is of the singlet type but, because of the term 
proportional to £±, its leading contribution is not only proportional to what one would expect 
from a pure singlet potential model. In Ref. [32] the authors have also elaborated on Eq. 
(152). Within the context of NRQCD [1], they use the leading equations of motion 9 , the 
power-counting rules of [1,25] and some arguments to neglect some mass-like terms, which 
could be generated under regularization. They get 

(^(^)|P 1 ( 3 5 1 )l^(^)>Ref. [32] = CA^^mE^l (157) 
(P Q (^)|P 1 ( 1 5 )|PQ(^)) Ref . [32] = C A ^f^mE^l (158) 



9 We have also used the equations of motion in order to derive Eq. (103). Nevertheless, we have 
done so in the context of pNRQCD. 
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where the term proportional to E\ is missing. Nevertheless, this does not necessarily reflect 
any inconsistency in any of the derivations since, according to the (perturbative-like) power- 
counting rules of [1,25], the term due to S± would be subleading. In any case, it would be 
very interesting to see how a term proportional to E\ would appear in the derivation of Ref. 
[32]. Here, we only would like to point out the possibility that an E\jm term may show up 
as a correction to the neglected mass-like term in Ref. [32]. Finally, let us note that in the 
dynamical situation mv ~ Aq CD , where £ 1 ~ Aq CD ~ m 2 v 2 ~ mE^d , both terms on the 
right-hand side of Eq. (152) are of the same order and contribute to the decay width at order 
c(a s (m))mv 5 . Phenomenologically this is particularly relevant to the case of pseudoscalar 
decays into light hadrons and to the electromagnetic decays. In the case of vector decays 
into light hadrons the contribution coming from the operator {VQ(nS)\Vi( 3 Si)\Vq(tiS)) may 
not be so important since the matching coefficient lmgi( 3 Si) ~ a s (m) 3 is suppressed by a 
factor a s (m) with respect to the other ones (with the exception of Im/ 1 ( 3 S' 1 ) and Im/ 8 ( 3 Pi), 
which are also of order a s (m) 3 ). 



D. Evolution equations 

In [1] evolution equations for the 4-fermion operators have been obtained. If we focus 
on the states that we are studying in this paper, the following evolution equations for the 
NRQCD matrix elements are obtained: 

(V Q (nS)\ LA-O^sS) \V Q (nS)) = ((V Q (nS)\O 8 ( 3 P )\V Q (nS)) 



dv ) 37vm 2 

+ (VQ(n5)|0 8 ( 3 P 1 )|VQ(n5)) 
+ (V Q (nS)\0 8 ( 3 P2)\V Q (nS)) 
-CfiVQinS^eS^VQinS))) , (159) 

(P Q (nS)\ [ u-fo^So)) \P Q (nS)} = -**L {(P Q (nS)\O s CPi)\P Q (nS)) 



dv J 37rm 2 

-CfiPQinS^CS^PQinS)}) , (160) 

(V Q (nS)\ (v^O^S!)) \V Q {nS)) = - 8 ^(V Q (nS)\V EM ( 3 Si)\V Q (nS)), (161) 

(P Q (nS)\ (^Oem(^o)) \P Q (nS)) = -^^(P Q (nS)\V EM CSo)\P Q (nS)). (162) 
Since we have, at 0(a s ) and leading log accuracy, 

u^-S 3 = 12C f — , (163) 

dv 7T 

Eqs. (147)-(150) are compatible with the evolution equations (159)-(162) at leading log 
accuracy. Note that at this order there is no v dependence in the states, and hence the 
derivatives with respect to v can be taken out of the expectation values. In Ref. [16] it was 
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proved that Eq. (163) gives the correct running for the octet operator of Eq. (156). In 
Appendix C, the reader can find the evolution equations and their leading-order solutions 
for the imaginary parts of all the 4-fermion matching coefficients needed in this work. 

VII. MODEL-INDEPENDENT PREDICTIONS 

The inclusive decays of the heavy quarkonium (either hadronic or electromagnetic) are 
usually considered up to, and including, NRQCD matrix elements of 4-fermion operators 
of dimension 8. This means to consider the (9(l/m 2 , 1/m 4 ) local 4-fermion operators of 
the NRQCD Lagrangian. With this accuracy, the decay into light hadrons of a vector S- 
wave state is described in NRQCD by the matrix elements of two singlet operators (0i( 3 Si) 
and V\), and three octet operators (0 8 ( 3 Si), O 8 ( 1 S , ) and 8 (P)). The corresponding pseu- 
doscalar S-wave state decay needs, at the same level of accuracy, the additional knowledge of 
the matrix element of the singlet operator Oi( 1 S'o). The electromagnetic decays of the same 
S-states need the additional knowledge of the matrix elements of the singlet electromagnetic 
operators Oem( 3 Si) and Oem( 1 *S'o) respectively. The decay of a P-wave quarkonium state 
into light hadrons and the corresponding electromagnetic decay are described in NRQCD 
with the above accuracy by the matrix element of a singlet (Oi(P)) and an octet (O 8 ( 1 5'o)) 
operator. If we consider that in the bottomonium system in principle 14 S- and P-wave 
states lie below threshold (T(nS) and i]b(nS) with n = 1,2,3; hb(nP) and Xb.j{nP) with 
n — 1,2 and J = 0, 1,2) and that in the charmonium system this is the case for 8 states 
(ip(nS) and i] c (nS) with n — 1,2; h c (lP) and % c j(lP) with J = 0, 1, 2), all the bottomonium 
and charmonium S- and P-wave decays into light hadrons and into photons or e + e~ are de- 
scribed in NRQCD up to (9(l/m 4 ) by 46 unknown NRQCD matrix elements (40 for the 
S'-wave decays and 6 for the P-wave decays) . These matrix elements have to be fixed either 
by lattice simulations [33] or by fitting the data [34]. Only in the specific case of matrix 
elements of singlet operators does NRQCD allow an interpretation in terms of quarkonium 
wave functions and one can resort to potential models. 

At the same level of accuracy S- and P-wave bottomonium and charmonium decays are 
described in pNRQCD, under the dynamical assumption Aqcd ^ rnv 2 , by only 19 non- 
perturbative parameters. These are the 13 wave functions (one for each of the 10 S'-wave 
quarkonium states below threshold, for which we need to distinguish different spin states, 
and a total number of 3 for the P-wave quarkonium states) and 6 universal non-perturbative 

parameters, which do not depend on the flavour and on the state (Si, S 3 , B\, S^ 2 '^ and 
£p' EM) ) 

In the above discussion we have counted NRQCD matrix elements by their dimensionality 
only. A more refined discussion would require that a maybe less conservative power counting 
be assigned to the NRQCD matrix elements as well as that the a 8 (m) suppression due to 
the short-distance NRQCD matching coefficients be taken into account. As we have already 
mentioned throughout the paper, the power counting of the NRQCD matrix elements is an 
open issue. To consider all the possibilities and phenomenological consequences goes beyond 
the scope of the present paper, whose aim is to set the theoretical framework. However, we 
would like to mention a few things. In the standard NRQCD power counting [25], the octet 
matrix elements are 0(v 4 ) suppressed for S-wave decays if compared with the leading order. 
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This is not so within our framework where, assuming the counting Aqcd ~ Tnv, they would 
only be 0(v 2 )-suppressed. This is potentially relevant to r(V — > LH) since Im/ 1 ( 3 5'i) is 
(9(o; s (m))-suppressed with respect Im fs(S). In other words, the octet matrix element effects 
could potentially be much more important than usually thought for these decays. It would 
be interesting to analyse this possibility further. 

The dramatic reduction in the number of parameters makes it possible, in the framework 
of pNRQCD, to formulate several new predictions with respect to NRQCD. In particular it is 
possible to relate information gained from decay widths of quarkonium with a specific flavour 
and principal quantum number to decay widths of quarkonium with different flavour and/or 
principal quantum number. Following this strategy in [16] the non-perturbative parameter 
£ 3 has been fixed on the charmonium P-wave decay data and used to predict ratios of P- 
wave decay widths for the bottomonium system (in this case and at leading order there is 
no ambiguity on the relative size between the singlet and the octet contribution). Here we 
will concentrate on some exact model-independent relations valid for S-wave decays. 

Let us consider the ratios of hadronic and electromagnetic decay widths for states with 
the same principal quantum number: 

_ T(V Q (nS) -+ LH) 
R "-T{V Q {nS)^e+e-y (165) 
P _ V{P Q {nS) -> LH) 

Ten of these ratios exist, ten being the number of bottomonium and charmonium states below 
threshold. As we discussed above, in NRQCD, and if one includes all the NRQCD operators 
up to (9(l/m 4 ), these 10 ratios are described by 40 non-perturbative parameters. It is a 
specific prediction of pNRQCD that, for the states for which the assumption Aqcd ^ mv 2 
holds, the wave-function dependence drops out from the right-hand side of Eqs. (165) and 
(166). The residual flavour dependence is encoded in the powers of 1/m, in E^j and in 
the Wilson coefficients, while the residual dependence on the principal quantum number is 
encoded in the leading order binding energy E^o ■ In principle, if all the 10 bottomonium and 
charmonium S-wave states below threshold belonged to the dynamical regime Aqcd 3> mv 2 , 
then, in the framework of pNRQCD, the ratios of hadronic and electromagnetic decay widths 
would be described by the 6 non-perturbative universal parameters listed above only. 

Particularly simple is, in pNRQCD, the expression of the ratios between R% and R% 
with different principal quantum number. We obtain up to order v 2 (with the counting 
Aqcd ~ mv) (M(nS) -2m = E$(l + 0(v)), M(nS) being the meson mass): 

Rl ( Img^S,) I^ \ IKH^ 

Rl + VWi( 3 Si) WeeC^ m 1 } 

K =l+ f lm<7i(%) _ Img^Sp) ^ M(nS) - M(mS) 

R p m VWi(%) lmf^S )J m 1 ' 

It is to be stressed that the octet-type contributions cancel (otherwise they would be l/a s (m) 
enhanced in the vector case). This prediction should be compared with the one expected 
in NRQCD. Within the standard (perturbative-like) power counting, the same prediction is 
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obtained in NRQCD. However, if one counts a s (m) ~ v 2 as was done in [32], the contribution 
due to the octet matrix elements is of the same order as the corrections obtained above and it 
should be taken into account in the vector case. Therefore, in principle, one is able to check 
the theory and/or the power counting. As an example, taking mi, = 5 GeV we get for the 
T(2S) and T(3S) states of the bottomonium system, Rj / Rj ~ 1.3, which is close (within a 
10% accuracy) to the experimental central value of about 1.4 that one can get from [35]. Let 
us also notice that, since Img 1 ( 1 S , )/Im i /' 1 ( 1 ,S , ) — Im g^i 1 S ) /lm f^i 1 S ) = 0(a 8 (m)), up to 
corrections of order v 3 we find that i.e. the ratio between hadronic and electromagnetic 
decay widths for pseudoscalar quarkonium, is the same for all radial excitations. However, it 
is not the purpose of this work to carry out a comprehensive and detailed phenomeno logical 
analysis, which is left to a subsequent publication. 

VIII. CONCLUSIONS 

We have obtained the imaginary part of the pNRQCD Hamiltonian up to 0(l/m 4 ) in the 
non-perturbative regime (k > Aq CD ^> mv 2 ). The expressions are given in Eqs. (87)-(89). 
As for any quantum-mechanical Hamiltonian, also the pNRQCD Hamiltonian is defined up 
to a unitary transformation. An alternative expression, related to the previous one by a 
unitary transformation, can be found in Sec. VD. 

We have applied our results to calculate the inclusive decay widths to light hadrons, 
photons and leptons up to 0(c(a s (m))mv 3 x (AQ CD /m 2 , E/m)) for S'-wave heavy quarkonium 
and up to 0(c(a s (m))mv 5 ) for P-wave heavy quarkonium. These are given in Eqs. (141)- 
(146) and are the main result of the paper. An alternative way to present it is given in 
Sec. VI C, where all the NRQCD matrix elements entering in quarkonium decays up to 
this order are expressed in terms of the quarkonium wave functions at the origin and 6 non- 
perturbative gluonic correlators, which are flavour- and state- independent, and for this reason 
may be called universal. The wave-function dependence factorizes in all these expressions. 
It is particularly remarkable that this is also true for the octet matrix elements. 

We have derived our expressions in two different ways: in Sec. Ill under the general 
assumption Aqcd < k and in Sec. V under the particular assumption k ^> Aqcd- In the 
first case, we have matched directly NRQCD to pNRQCD in an entirely non-perturbative 
one-step procedure, based on the Hamiltonian formulation of NRQCD. In the second case, 
we have matched NRQCD to pNRQCD in a two-step procedure, the first perturbative, 
the second non-perturbative, but still with a clear diagrammatic interpretation based on 
the multipole expansion. The fact that these two completely different ways of deriving 
the pNRQCD Hamiltonian give the same answer up to a unitary transformation can be 
considered a stringent test on the correctness of the result. In Sec. VI D we have also 
checked that the evolution equations of our universal parameters are consistent at leading 
log accuracy with the known evolution equations of the NRQCD matrix elements. 

In Sec. VII we have considered the phenomenological implications of our results. There 
exist 14 charmonium and bottomonium states below threshold. We expect our results to be 
applicable to most of these states. The exceptions are, on the one hand, the T(IS'), which is 
commonly understood as a weak coupling state (i.e. k ^> E > Aq CD ), and, on the other hand, 
states that are too close to the D-D threshold for charmonium or to the B-B threshold for 
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bottomonium, like, maybe, the ip(2S). Going from NRQCD to pNRQCD reduces the number 
of non-perturbative parameters needed to calculate the inclusive decay widths associated 
with these states by about a factor of 2. The situation is even better if we consider ratios of 
hadronic and electromagnetic decay widths. Since the wave-function dependence factorizes, 
it drops out in the ratios. It follows that only 6 universal parameters, which depend only on 
the light degrees of freedom of QCD, are needed. The already known data will be sufficient 
to fix all these parameters, to allow checks and to make new predictions. Moreover, suitable 
combinations of ratios give rise to novel parameter- free, model-independent predictions. We 
have considered some of them in Sec. VII. 

The non-perturbative universal parameters that we have introduced do not necessarily 
need to be fitted to the experimental data. We have provided expressions for them in 
terms of correlators of gluonic fields. This allows for an eventual evaluation on the lattice. 
These parameters may also be obtained from QCD vacuum models [36]. We note that, 
once they become fixed, our results make the evaluation of NRQCD octet matrix elements 
possible from properties of the wave functions at the origin. Hence, any potential model that 
leads to definite wave functions [37] will provide definite results for these matrix elements. 
Nevertheless, it should be pointed out that, if we wish to obtain the NRQCD matrix elements 
given in Eqs. (147)-(150) with the aforementioned accuracy, any potential model to be used 
here must be consistent with the structure of the potential derived from NRQCD in terms 
of Wilson loops in Refs. [14,15]. In fact, the wave functions defined in this paper can 
also be computed in a model-independent way without resorting to data fitting. This is 
so because our wave functions correspond to the solution of a Schrodinger equation where 
the potentials are given in terms of expectation values of Wilson loops with suitable field 
insertions. Therefore, once lattice simulations are provided for the potentials [38], the wave 
function can be obtained unambiguously without any model dependence. 

Since our method reduces the number of unknown parameters with respect to NRQCD, 
we expect it to become increasingly relevant as the number of needed NRQCD matrix ele- 
ments increases. This seems to be necessary in the calculation of charmonium decay widths, 
where the non-relativistic expansion converges slowly. Indeed, higher order operators have 
been considered recently in Refs. [39,40]. In Appendix B, we give the general matching 
formula for the NRQCD matrix elements to the pNRQCD results without going through the 
whole matching procedure outlined in the main body of the paper. 

We have also addressed, mainly in Sec. IV B, the issue of the power counting in NRQCD 
in the non-perturbative case. We believe that our formalism provides a suitable theoretical 
framework to study it. The power counting of NRQCD is not known a priori in the non- 
perturbative regime and it could, in principle, be different, depending on each dynamical 
system. This is particularly transparent in pNRQCD. There, the potentials are functions of 
r and Aqcd- Therefore, as the typical value of r changes from system to system, one should 
accordingly assign a different size to each given potential. Moreover, having expressed the 
NRQCD matrix elements in terms of wave functions and universal correlators, we have 
disentangled the soft scale k, now entering in the wave function square, from the Aq CD /m 
and E/m corrections. In fact, this is why we can construct ratios of convenient decay rates 
where the /c-dependence drops, providing a more constrained set of relations. For these ratios 
the fixing of the power counting reduces to the evaluation of the correlators, while taking 
into account possible enhancement effects due to the NRQCD matching coefficients. 
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Finally, although in the present paper we have focused on inclusive decays to light 
hadrons, there should be no conceptual problem, a priori, in considering the NRQCD matrix 
elements that appear in heavy quarkonium production. We also expect there a significant 
reduction in the number of non-perturbative parameters. In particular, our formalism may 
shed some light on the power counting problems that appear in the heavy quarkonium po- 
larization data [28]. 
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APPENDIX A: 4-FERMION OPERATORS 



[1], 



Here we list the relevant 4-fermion operators of dimension 6 and 8, as taken from Ref. 
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t/ 



i/; T 
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2 L 
1 

2 
1 
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-D x a)x-xK-0 x 
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',t v v t/ 



^ T xX T (-fD)> + H.c. 
^o-x-xM-lD) 2 ^ + H.c. 
V>VxxV(-i) 2 D (i iyty + H.c. 



V — 2^;^ A A V 2 / ^' 



- V f (-|D x o-)T a X • x f (-|D x cr)T a ip, 



, 2 ^ " 7 AAV 2 ^ " 7 J 

^XX^-fD) 2 ^ + H.c. 
^ f oT a x • x t o-(-|D) 2 T a -0 + H.c. 
V>VT a xx t ^'(-f) 2 D (i D J ')T a -0 + H.c. 



(Al) 
(A2) 
(A3) 
(A4) 

(A5) 
(A6) 

(A7) 
(A8) 
(A9) 

(A10) 

(AH) 
(A12) 
(A13) 

(A14) 
(A15) 
(A16) 

(A17) 

(A18) 



where we use the conventional notation T fe) = (T ij + T ji ) /2-T kk 5 ij /3. The electromagnetic 
operators are defined as follows: 



OemCSo) 


= ^ f x|vac) (vac|xV, 


(A19) 


Oem( 3 S\) 


= ^Vxlvac) (vac|xW, 


(A20) 


OemCPi) 


= ^(_iD) x |vac) • (vac| X t(_iD)^, 


(A21) 


Oem( 3 Po) 


= 1 V f (-|D • <T)x|vac>(vac| X t(-iD • tr)^, 


(A22) 


O EM ( 3 Pi) 


= \ ^(-|D x <r)x|vac) • (vac|x f (-|D x <r)^, 


(A23) 
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E m( 3 P 2 ) = V>V|D ( V))x|vac)(vac|x^fD(V'))V>, 



VemCSq) — - 



^em( 3 5'i] 



■0^<rx|vac) (vac^cr ( — -D 2 ) ■?/> + H.c. 



where |vac) is the vacuum state of QCD. 



(A24) 
(A25) 

(A26) 
(A27) 



APPENDIX B: DIRECT MATCHING TO PNRQCD OF NRQCD MATRIX 

ELEMENTS 

In principle, it is possible to match directly to pNRQCD matrix elements of NRQCD that 
involve operators different from the Hamiltonian H. In this way NRQCD matrix elements 
can be expressed in terms of non-local correlators without going through the full matching 
procedure outlined in the main body of the paper. This is useful if no iteration of these 
NRQCD operators are necessary in the matching calculation. In order to do this it is 
necessary to have an explicit expression for the state |0;xi,x 2 ). Up to 0(l/m) it can be 
found in Eq. (36). This way of doing will be particularly useful in order to work out higher 
order operators that will appear in going beyond 0{mv b ) in the expansion of the heavy 
quarkonium decay width. Higher order operators appear to be necessary for charmonium 
decays, where the non-relativistic expansion converges slowly, assuming v 2 c ~ 0.3. 

The master equation is (\H) represents a generic heavy-quarkonium state at rest, P = 0, 
with quantum numbers n, j, / and s as defined in Ref. [1]): 



(H\0\H) 



(P = 0|P = 0) 



J d 3 r J d 3 r' J d 3 R J d 3 R' (P = 0\R){njls\r) (Bl) 



(0; Xl x 2 | d 3 ZO(Z) |0;x' lX2 > (R'|P = 0){r'\njls) , 



where r = xi — x 2 , r' = 
(R'|P = 0) = 1 and (P 
NRQCD matrix element 



x 



< - x' 2 , R = (xi + x 2 )/2 and R = (xi + x' 2 )/2 (note that 
0|P = 0) = J d 3 x). As an example, let us consider here the 



( XQ (n01)\F EM ( 3 Po)\XQ(n01)), 



of the dimension-9 operator 



Fem( 3 Po) = ~^<T ■ ^ExIvacXvaclxV • + H.c, 



(B2) 



(B3) 



which is relevant to describe the electromagnetic decay Xco 77 at order mv 7 accuracy. 
Owing to spin symmetry, the same matrix element enters into the Xc2 77 decay. These 
contributions have recently been considered in [39]. In the Hamiltonian formalism of Sec. 
Ill the matrix element (B2) is written as 
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(xo("Oi)iy EM ( 3 -Po)ixo("Oi)) = - | p - 2 y y dvy y ^p^ir) 

x(„011|r) [y^fex.x /^-^^'^'^ ffl |0;x<^)" 

x (r'|n011)(R'|P = 0), (B4) 

where |n011) is the Schrodinger wave function of the state xq(^01). Now we expand the 
state (0;X]X2| according to Sec. Ill C. The first non- vanishing contribution comes from the 
1/m correction given in Eq. (36). Inserting that expression into Eq. (B4) and keeping in 
mind that only the term with the derivative projects onto the |n011) state, we obtain 

(XQ(n01)|^EM( 3 Po)|XQ(n01)) = Q | p - - 2 J d 3 r J d 3 r' J d 3 R J d 3 R' (P = 0|R) 

x(°)(k;x lX2 | I ^ ^■^x|vac)(vac| X t (T .D^ ^ ^ (0) ^ , = 



2 



E (»>(0| g g»'°'^(%g|va C > (va C |0)(°) ( „ 011|j , v ,,( 3)(r)(r . v|nOU) (B5) 



3m ^ (£<°>-£f>« 



In the second equality we have made use of Eq. (53) and of the Wick theorem. Finally from 
the fact that 5( 3 )(r)|0)W = 5( 3 )(r)l c /v^Vc|vac) and from Eq. (73) we get 

(XQ(nOl)\F EM ( 3 Po)\XQ(n01)) = -C A ^ nl Wl ^, (B6) 

it m 

or equivalently, using Eq. (151), 

(XQ(n01)|^ EM ( 3 Po)|XQ(n01)> 2S 1 



m( XQ (n01)\0 EM ( 3 Po)\XQ(n01)) 3 m 



2 ' 



(B7) 



Similar considerations may in principle also be applied to the matrix elements needed at 
relative order v 4 for S'-wave decays. For a complete set of these operators and for consider- 
ations concerning their relevance in phenomenological studies, see Ref. [40]. 



APPENDIX C: RUNNING EQUATIONS OF THE MATCHING COEFFICIENTS 

The running equations obtained in Appendix B.3 of Ref. [1] for the NRQCD 4-fermion 
operators give us information on the running of their matching coefficients. The running 
equations read as follows 

vflmnfSo) = ICf^lmtfSo) , (CI) 
av 6 7i 

is-^-lmg 1 ( 3 S 1 ) = ^C>— Im/i( 3 Si) , (C2) 
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and otherwise. 

The imaginary pieces of the dimension-6 operator matching coefficients {f{S)) do not 
run at leading non-vanishing order: 



^Im/(S) = 0. 



(C15) 



Therefore, the above equations can be easily solved in that case. We obtain at leading 
non-vanishing order: 
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where we have chosen to as the starting point of the evolution equation; the matching 
conditions at this scale at 0(o%) can be read from Ref. [1]. 

APPENDIX D: REGULARIZING PRODUCTS OF DISTRIBUTIONS 

In the intermediate steps of the calculation we find ill-defined products of distributions. 
We first show how dimensional regularization makes sense out of these expressions by setting 
them to zero, and next how they amount to renormalizations of local terms when a cut-off 
regularization is used instead. 

Consider, first, the product of two delta functions: 

5(3) ( r ) S<?\r)=J d D p J d D p' J rfV|p)(p|5 (3) (r)|p / )(p , |^(r)|p")(p ,/ | 

= J d D p j d D p' J d D p" |p)(p"| = 0, (Dl) 

since the integral over p' has no scale. 
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Consider next: 



* <3) ( r )^ = J d °P / dD P'J rf D P"lp)(pl« ,3, (r)|p'><p'l^|p"><p"l 

= Jd° P J d° P 'J d° P " Ip> |p ,^ |3 -, (P"I = 0, (D2) 

since, upon the translation p' — > p' + p", the integral over p' has no scales. 

Alternatively, if we use a cut-off regularization, for instance by smoothing the delta in 
momentum space, like: 



(p|^ 3 )(r)|p'> = 1 ^ e -^, (D3) 



we obtain: 



5^(r) 8^(v) ~ ^La^)^) - ({V 2 , ^(r)} + 2V^ 3 )(r)V i ) + O , (D4) 

which can be removed by local counterterms. Hence DR implements nothing but a suitable 
subtraction prescription. Analogously, it is easy to see that S^(r)/r s for s = 0,1,2,... 
reduces to local terms. 



APPENDIX E: UNITARY TRANSFORMATIONS 

It is well known that quantum-mechanical Hamiltonians, which are related by unitary 
transformations, lead to the same physics. This fact is particularly relevant to quantum- 
mechanical Hamiltonians, which are derived from a field theory, which is our case. It is also 
well known that the quantum-mechanical potentials that are obtained from QED depend 
on the gauge one uses in the calculation (this is also so for QCD in perturbation theory), 
but physical observables computed with either potential turn out to be the same. It is 
perhaps not so well known that the potentials obtained in one gauge can be related to the 
ones obtained with a different gauge by means of a unitary transformation. In fact the 
arbitrariness in the form of the potentials is not only due to gauge dependence. It depends 
in general on the way one carries out the matching calculation. Any correct result is related 
to any other one by means of a unitary transformation. 

We shall use this fact here to prove that the result obtained in Sec. V D is equivalent to 
the one obtained in Sec. IV A. 

Consider the following unitary transformation: 

U = e l{r ' p} ^ . (El) 
Consider next a delta function in the Hamiltonian: 

m 2 m 2 m 6 L J 

which on S-wave states reduces to: 
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This shows that a suitable unitary transformation may induce terms at 0{\/m 3 ) proportional 
to |0(O) | 2 . Of course, physics should not change. If 0(r) is an eigenfunction of h, then 
0(r) = W(j)(r) is an eigenfunction of h = WhU. Then: 

0(0) = e- i(2rp - 3t) ^0(O) (1 + 0(|r|)) | r=0 = e~ 3 ^0(O) . (E4) 

Clearly: 

^ + ^) \m\ 2 - — 2 \m\ 2 ■ ( E5 ) 

We will illustrate this issue further with an example. Recall the two different results, 
namely (116) and (117), we obtained from the first diagram of Fig. 2 at second order in the 
expansion ?to; 2 /Aqcd- More explicitly, from (116) we get a real result: 



4V Z 



E — h s m 2 

and from (117) a result containing an imaginary part: 



f°° 1 
J dtt'igE^-gEiO))^—^, (E6) 



-2(t.™+W(™ +1 Ss)^) rdtt 2 ( 9 E(t). 9 E(0))-^—. (E7) 
> \ m m 6 J J E - h s 



E-h s 

Both results are correct. They are indeed related by the following unitary transformation: 

U = e^'^ , 
1 f°° 

£2= N~cJo ^ ' ^ E(0)) ' (E8) 

and hence lead to the same decay width. This can even be further confirmed by an explicit 
calculation in the case of the Coulomb potential, since the induced terms then retain the 
same form as the original ones. 

The unitary transformation, which brings the result of Sec. IV A to the one of Sec. VD 
reads: 



U = e 



^ = + 4 2 > EM ) _ Eft) _ gp™)) . ( E9 ) 

Clearly this transformation also reshuffles 1/m real potentials into 1 / m 3 real potentials. 
This means that, in the more conservative counting considered in Sec. IV B, the whole set of 
potentials up to 0(l/m 3 ), which are formally given in Sec. Ill B, are expected to be relevant 
to calculate the wave function at the origin with an accuracy that matches the NNLO terms 
in (141)-(146). 
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FIGURES 




Figure 1 : The interaction vertices in pNRQCD' which are needed in order to calculate the decay width up to ^ . 
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a) 




b) 



Figure 2: Diagrams corresponding to the leading contributions in ( A ^° D ) (Fig. a)) and ( A< ^ £> ) 
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(Fig. b)). All corrections not contained in £ 3 
propagators. This generates all terms exhibited in Figure 3 



and £P arise from them after expanding the internal 
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Figure 3: Diagrams generated by those in Fig. 2. a) corresponds to a P wave octet correction; b) and c) give rise 
to a chromomagnetic two-field correlator accompanying both a spin-flip/octet and a non-flip/singlet imaginary 
coefficient, respectively; d) produces the term proportional to £ s times the binding energy; e) shows the structure 
introduced by the lmgx ( 2S+1 S s )-proportional contact interaction. 
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Figure 4: Diagrams corresponding to the next-to-leading contributions in the ^ j expansion. After expan- 
sion of the internal propagators, as explained in the text, they produce the series of graphs presented in figures 5 
to 9, which originate the terms proportional to £f'^ and £^ 2 ' . 
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Figure 5: Diagrams stemming from Fig. 4a) and 4b). They arise from terms of the form 
i , (h s — E)] | ([(h s - E), ] + [ ,(h s — E)])[(h s — E), ] upon expansion of the octet propagators. 
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Figure 6: Diagram generated by Fig. 4c) after picking up in the expansion of the middle octet propagator the 
termof the form [ ,(h s - E)] ((V (r) - V s (r)) [(h„ - E), }. 
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Figure 7: All remaining diagrams generated by Fig. 4c). Here the expansion of the octet propagators keeps 
the sequence [ , (h s - E)] | ({(h s - E), ] + [ , (h s — E)]) {(h s - E), ], with suitable gluonic vertices in- 
serted in each case. 
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Figure 8: Diagrams generated by Fig. 4d) after projecting out the vacuum insertion and upon expanding both 
singlet and octet propagators. As seen, gluonic vertices are conveniently inserted in a propagator sequence of the 
form[ ,{h a -E)]\{[{h a -E), ] + [ , (h s — E)]) [(h s — E), ]. 
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H:= -i{h s - E) 



Figure 9: Diagrams contributing to the potential generated by the vacuum insertion in Fig. 4d). a) causes the 
structure £480 to appear. The four of them are responsible for the £381 term. The operators A and ■ act through 
suitable commutations, which are not reflected in the figures, on the vertices. "H must be taken left and right 
according to the prescription given in the text. 
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